Gas Properties and Implications for Galactic Star Formation in Numerical 
Models of the Turbulent, Multiphase ISM 



Hiroshi Koyama 1,2 and Eve C. Ostriker 1 
ABSTRACT 
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Using numerical simulations of galactic disks that resolve scales from ~ 1 to several 

hundred pc, we investigate dynamical properties of the multiphase interstellar medium 

in which turbulence is driven by feedback from star formation. We focus on effects 

Q ' of HII regions by implementing a recipe for intense heating confined within dense, 

self-gravitating regions. Our models are two-dimensional, representing radial-vertical 

slices through the disk, and include sheared background rotation of the gas, vertical 

stratification, heating and cooling to yield temperatures T ~ 10 — 10 4 K, and conduction 
<**] , 

that resolves thermal instabilities on our numerical grid. Each simulation evolves to 
reach a quasi-steady state, for which we analyze the time-averaged properties of the 
gas. In our suite of models, three parameters (the gas surface density S, the stellar 
volume density p*, and the local angular rotation rate Q) are separately controlled 
in order to explore environmental dependencies. Among other statistical measures, we 
evaluate turbulent amplitudes, virial ratios, Toomre Q parameters including turbulence, 
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and the mass fractions at different densities. We find that the dense gas (n > 100 cm" 
has turbulence levels similar to those observed in giant molecular clouds and virial ratios 
~ 1 — 2. Our models show that the Toomre Q parameter in the dense gas evolves to 
values near unity; this demonstrates self-regulation via turbulent feedback. We also 
test how the star formation rate Ssfr depends on S, p*, and Under the assumption 
that the star formation rate is proportional to the amount of gas at densities above a 
threshold n t ^ divided by the free-fall time at that threshold, we find that Ssfr oc 
with 1 +p ~ 1.2 — 1.4 when nth = 10 2 or 10 3 cm -3 , consistent with observed Kennicutt- 
Schmidt relations. Estimates of star formation rates based on large-scale properties (the 
orbital time, the Jeans time, or the free-fall time at the mean density within a scale 
height), however, depart from rates computed using the measured amount of dense 
gas, indicating that resolving the ISM structure in galactic disks is necessary to obtain 
accurate predictions of the star formation rate. 

Subject headings: galaxies: ISM — hydrodynamics — ISM: general — method: numer- 
ical — instabilities, turbulence — stars: formation 
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Introduction 



The interstellar medium (ISM) is commonly envisioned as a self-regulating system in statisti- 
cal quasi-equilibrium. Multiple components of gas with varying densities and temperatures coexist 



(Field et al 



1969 



Cox Smithlll974l; iMcKee &: Ostrikerlll977l ). animated by turbulence that per 



vades the whole volume (jElmegreen &: Scaldl2004l ) . Different components of gas play different roles 
in the ISM ecosystem, with the coldest and densest portions responsible for star formation. Mas- 
sive stars, when they are born, energize the ISM through the HII regions and supernova blasts 
they create (|Spitzerlll978l ): this energy input is important in replenishing continual losses through 
turbulent dissipation. UV radiation from young massive stars is also crucial in heating the gas. 
The rate of star formation is determined by the available supply of dense gas, which in turn is reg- 
ulated by the interplay between dynamics and ther modynamics in the ISM, and is affected by the 
galac tic environment in which the ISM is contained (|Mac Low &: Klessenl 12004 ; IMcKee &; Ostriker 
20071 ). While this overall framework is generally accepted and is supported by existing theory and 
observations, much work remains on both fronts to quantify the dependence of statistical properties 
on the global system parameters, and to establish when and how self-regulated quasi-steady states 
are achieved. 

Given the importance of time-dependent processes and interdependencies in the ISM, complex 
theoretical models are needed in order to address even rather basic questions. For example, what 
sets the r elative proportion s of the different gas components? In an idealized classical picture such 
as that of lField et al.l (|1969l ) for the atomic medium, given a pressure and a mean density n, thermal 
equilibrium defines a density for each of two stable phases, n warm and n co \a, and the ratio of cold 
to warm gas is given by a simple algebraic relation: M co id/M warm = (n~ arm — n _1 )/(n" 



n 



-1 



cold J 



In the real ISM, however, which is a time-dependent system, thermal equilibrium only holds to 
the extent that the radiative times are short compared to dynamical times for compressions and 
rarefactions. Furthermore, the value of the mean density n and pressure (averaged over large 
scales) are not even known a priori for a given ISM surface density, since the vertical distribution 
of gas is sensitive to its dynamical state. This dynamical state itself depends on the (unknown) 
dense gas fraction, since more dense gas produces more feedback from star formation, and hence 
more turbulence that inflates the disk vertically to reduce n (and also produces local variations in 
density and pressure through compressions and rarefactions). Multidimensional effects (the ISM 
is not simply stratified perpendicular to the galactic plane, but is composed of filamentary clouds) 
and self-gravity additionally complicate the situation. 

In recent years, a number of groups have begun developing models of the turbulent, multiphase 



star formation (e.g.. Korpi et al.lll999: de Avillez & Breitschwerdt 200^ 


. 2005 


3. 2007 


Mac Low et al. 


2005: Joung & M 


ac Low! 2006. 


Dib et al. 2006). self-gravity of the gas (e.g.. 


Wada & Normal 1999: 


Wada et al. 2002 


). and both of these effects (e.g.. Wada & Norman 


2001, 


2007; 


31vz et al. 


2005; 


Tasker & Brvan 


2006. 2008: ] 


Robertson & Kravtsov 2008). The treatment of feedback in these 



simulations is to inject thermal energy in regions identified as sites of star formation; most models 
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focus on the energy input from supernovae. In very large scale simulations that have minimum 
resolution of only 50-100 pc, feedback implemented via thermal energy deposition is not in practice 
very effective, because the input energy is easily radiated away. With finer numerical resolution, 
feedback regions expand adiabatically at first to make hot diffuse bubbles, driving shocks that sweep 
up surrounding gas and ultimately generate turbulence throughout the computational domain. A 
number of different issues have been addressed by these recent simulations, including investigating 
departures from thermal equilibrium in density and temperature PDFs, measuring the relative 
velocity dispersions of various gas components, and testing whether relationships between star 
formation and gas surface density emerge that are similar to empirical Kennicutt- Schmidt laws. 

Even before the advent of supernovae, massive stars photoionize their surroundings, creating 
HII regions within molecular clouds that are highly overpressured and expand. HII regions may 
in fact be the most important dynamical agents affecting the properties of dense gas in giant 
molecular clouds (GMCs), since the original turbulence i nherited from the diffuse ISM is believe d 



to dissipate within a flow crossing time over the cloud (IStone et al 



1998; 



Mac Low et al.lll998h 



while GMCs are thought to live for at least a few crossing times (IBlitz et al.l 120071 ). Analytic and 
semi-analytic treatments find that GMCs with realistic sizes, masses, and star formation rates can 
indeed be maintained by the energy input from HII regions for a few crossing times, ultimately 
being destroyed thro ugh a combination of photoevaporat i on and kinetic energy inputs that unbind 



the remaining mass (IWhitworth 



1979 



Franco et al 



1994 



Williams McKee 



1997; 



Matznerll2002l 



Krumholz et al.l 120061 ) . Recent three-dimensional numerical studies have begun to ad dress this 
process in detail (e.g. iMellema et al.l 120061 : iMac Low et al.1 120071 : Krumholz et al.ll2007l ). focusing 
on regions within GMCs. 

In the present work, we consider how the large-scale dynamical state of the ISM is affected by 
star formation feedback in the form of expanding HII regions. Our main interests are in exploring 
how the turbulence driven by HII regions affects the properties of dense gas (we measure statis- 
tics of density, temperature, and velocity), in testing ideas of global self-regulation by feedback 
(we evaluate Toomre Q parameters and virial ratios), and in exploring how galactic environment 
systematically affects the character of the ISM, including its ability to form stars. Complete ISM 
models should of course include feedback from supernovae as well as those from HII regions, and 
it is our intention to do this in future work. However, we consider it useful to adopt a sequen- 
tial approach, independently testing the effects of HII region feedback to provide a baseline for 
more comprehensive simulations. In addition to developing a physical understanding of the ways 
in which feedback affects the ISM, another goal of our work is to investigate the sensitivity of 
numerical results to prescriptions that are a necessary - but not always fully tested - aspect of 
galactic-scale studies of star formation. In particular, we examine how the choice of density thresh- 
old in commonly-adopted recipes for star formation affects the resulting dependence of the star 
formation rate on ISM surface density. 



Our approach to exploring the effects of galactic environment is to conduct a large suite of 
local simulations that cover a range of values for three basic parameters: the total surface density 
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of gas in the disk (£), the local midplane stellar density (p*), and the local rate of galactic rotation 
(ft). The parameter range covers a factor of six in gas surface density and galactic angular rotation 
rate, and a factor of 30 in stellar density. Our suite is divided into four series, each of which has one 
independent parameter that is systematically varied. We also include comparisons with hydrostatic 
models that are identical in terms of their input parameters to the fully-dynamic models, but do 
not include feedback and hence are not turbulent. For this first set of pilot studies, we have not 
implemented full radiative transfer to evaluate the extent of HII regions (we intend to do so in 
the future) , but instead introduce a simple prescription in which the boundaries of HII regions are 
determined by the gravitational potential. Using this approach (rather than, for example, adopting 
a single fixed outer radius) has the advantage that the volume of the heated region expands as 
the density surrounding the source drops. Since our treatment of HII regions does not attempt 
to be exact, we do not consider our specific results for e.g. velocity dispersions to be more than 
approximate (although in fact we find similar values for velocity dispersions in dense gas to those 
that are observed in GMCs). Instead, we shall emphasize the general properties of a multiphase 
ISM system in which turbulence is driven from within the dense phase. 

This paper is organized as follows: In §2 we describe our numerical methods, and in particular 
the recipe for star formation feedback. The control parameters for our disk models, and the 
properties of each model series, are presented in §3. Section 4 gives an overview of evolution based 
on our fiducial model. In §5 we present the statistical properties of the gas in each model, and test 
environmental influences by intercomparing the model series. The implications of our results for 
star formation, both in real galaxies and in numerical simulations, is analyzed in §6. We conclude 
with a summary and discussion in £J7J 



2. Numerical Methods 



2.1. Basic Equations 



We study the evolution of rotating, self- gravitating, galactic gas disks, including local heating 
and cooling terms. We solve the hydrodynamic equations in a local Cartesian reference frame 
whose center lies at a galactocentric radius Rq and orbits the galaxy with a fixed angular velocity 
fto = ft(-Ro)- In this local frame, radial, azimuthal, and vertical coordinates are represented by 
x = R—Rft , y = Rp((j)—tlot), and z, respectively, and terms associate d with coordinate curvature are 
neglected (IGoldreich &: Lynden-Belll Il965bl ; lJulian &: Toomrd Il966l ) . The local- frame equilibrium 
background velocity relative to the center of the box at x = y = z = is given by vo = —qQoxy, 
where 



(fin ft 
dlnR 



(1) 



ft 
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is the local dimensionless shear rate. In terms of q, the local epicyclic frequency k is given by 

K2 ^^jR {R4n2) = {A - 2q)n2 - (2) 

We shall choose q = 1 for all models, representing a flat background rotation curve V c = £IR = 
const, for the unperturbed motion. 

In addition to the tidal gravity and Coriolis terms from the "shearing sheet" local formulation, 
we also include terms for the vertical gravity of the stellar disk, gas self-gra vity, radiative heat 



ing and cooling, and the rmal conduction. The resulting equations (see e.g., lHawlev et al 



1995 



Piontek fc Ostrikeri 12004 ) are: 



|f + V-(pv) = 0, (3) 
dv 1 

— + v-Vv = — VP + 2gOarx - 2fi X v- V$ + g*. (4) 
at p 

— + V • (ev) = -PV • v + n [T(t, x) - nA(T)] + KV 2 T, (5) 

where P = fnk^T and n = p/u. With n the number of hydrogen nuclei per unit volume, / varies 
from 0.6 to 1.1 depending on whether the gas is predominantly molecular or atomic; we simply 
adopt / = 1.1. We adopt u = lAm p to include the contribution of helium to the mass density. 
Here e = P/(7 — 1) is the internal energy per unit volume (we adopt 7 = 5/3), K is the thermal 
conductivity, $ is the self-gravitational potential due to gas, and the vertical gravitational force 
due to stellar disk is 



g* = —4irGp*zz, (6) 

where p* is the stellar density and z is the vertical coordinate relative to the midplane. In the 
above expressions and elsewhere in the remainder of the text, we have dropped the "0" subscript 
on k also refers to the value evaluated at the center of the domain. Computation of the gas 
self-gravity is discussed below. 

In this paper, we present results of two-dimensional simulations of the above set of equations. 
The two independent spatial coordinates in our models are x and z; thus, we follow evolution 
in radial-vertical slices through a galactic disk. Although the azimuthal (y) direction is not an 
independent spatial variable for the current set of models, we do include azimuthal velocities, and 
their variation with x and z. Inclusion of v y is important because angular momentum can strongly 
affect the ability of self-gravitating perturbations to grow. Radial motions that are required for gas 
to become concentrated are coupled to azimuthal motions through the Coriolis force; perturbations 
in the azimuthal velocity with respect to the mean background shear correspondingly lead to radial 
motions via the Coriolis force. Although our two-dimensional models do capture some of the 
effects of galactic rotation (i.e. epicyclic oscillations), they miss some of the effects associated 
with shear. In three dimensions (or in the height-integrated R — <j> plane), azimuthal shear can 
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make it more difficult for self-gravitating concentrations to grow. Of course, in three dimensions, 
self-gravity also increases more rapidly as the density increases, which enhances the ability of 
dense concentrations to grow. Although it will be important to revisit the present models with 
fully three-dimensional simulations, we do not anticipate large changes based on dimensionality. 
Previous three-dimensional simulations of shearing, rotating disks have found similar (within a 
factor 2) nonlinea r instability thre s holds for self-gravitating clou d formation to reduced-dimensional 
models (see e.g. iKim et al.l I20021 . 120031 : iKim fe Ostrikerl 120071 . and references therein). Thermal 
instability also develops similarly in the t wo-dimensional and three-dimensional case to crea te a 
cold cloud/warm intercloud structure (e.g. iPiontek fc Ostrikerl 120041 ; iPiontek fc Ostrikerl 12005). 



2.2. Hydrodynamic Code and Boundary Conditions 

The numerical solutions to the two-dimensional dynamical equations are obtained using a 
temporally and spat ially second order finite volume method which includes TVD Runge-Kutta 
integration in time (jShu &; Osherl Il988l ). with a directional ly unsplit flux update and piecewise 
linear reconstruction with slope limiter (see, e.g.. lHirschll2007l ). We use Roe's approximate Riemann 
solver with an entropy fix (jRogJ |l98l|) . The heating and cooling terms in the energy equation are 
separated out in an operator-split fashion and updated using implicit time integration (see §2.4j) . 
because the cooling times are frequently much shorter than the other timescales. The code is 
parallelized using MPI. 

For the hydrodynamic update, the time step is set to At = min(iHD> ^cond) *cooh *heat) where 



^HD 

^cond 
^cool 
^heat 



Ccfl min 

Ck min 
Ct min 
Ct min 



1 



n/c B (Ax) 2 



( 7 -l)nA(T) /all zoncs 



( 7 -i)r(x)y alI 



(7) 

(8) 
(9) 
(10) 



and we adopt Ccfl = 0.8, Ck = 0.5 and Ct = 50. Here, c s = (jkT/ fi) 1 / 2 is the sound speed in 
any zone. With a large value of Ct, the explicit hydrodynamic timestep is not strongly limited by 
the cooling time in dense gas. The adopted Ct is chosen such that the solution agrees with tests 
of individual expanding "HII regions" (for our feedback model) that have Ct ~ 1 (equivalent to 
explicit cooling); if a much larger value of Ct is allowed, this expansion is not accurately reproduced. 



At the x (radial) boundaries, we implement shearing-periodic boundary conditions (jHawlev et al 
19951 ). in which the azimuthal (angular) velocity term v y = R${4> — Oo) = (Ro/fyv^ — V c is in- 
cremented or decremented by ±L x £Iq in mapping from the right— ► left or left— > right boundary, 



-7- 



consistent with the equilibrium velocity field. In the z-direction, we adopt periodic boundary con- 
ditions for the hydrodynamic variables, such that the total mass in the domain is conserved^ The 
gravitational potential solver applies open (i.e. vacuum) boundary conditions in z, as we next 
discuss. 



2.3. Poisson Solver 

We have developed a new method for obtaining the gravitational potential of a disk in Cartesian 
geometry using Fast Fourier Transforms (FFTs). Since the discrete Fourier Transformation allows 
only periodic functions, a special approach is needed to solve for a disk potential with vacuum 
boundary conditions outside the simulation domain. 

Let us consider a simple case, consisting of an a uniform, isolated gas sheet in the z = plane 
which has density p(z) = a5{z). The corresponding gravitational potential <&(z) = 2nGa\z\ is 
obtained by solving the Poisson equation, 

V 2 $ = 4vrGp, (11) 

with vacuum boundary conditions. If we have a finite domain of size L z and suppose that the 
gas sheet lies somewhere within the domain, then we only would require values of the potential at 
locations within ±L Z of the sheet. Thus, we may find the potential within z = (—L z , L z ) in terms 
of discrete Fourier components as 

$ e = 2irGaC t , (12) 

c = [ L - ]zl M^ = - 2 J^p. d3) 

J-l z (si 

\L-, 

In Fourier space, the Poisson equation for one independent variable is 

- k 2 $> k = 47rGp k . (14) 

Thus, in terms of discrete Fourier components with hi = £(2ir /2L Z ), we have 

i 4xGpt . . 

®e = —, — 72- ( 15 ) 



We have found that except for mass loss, the overall evolution is similar when we apply outflow boundary 
conditions in the vertical direction. Adopting periodic boundary conditions for hydrodynamic variables makes it 
possible to maintain the gas surface density E at a constant value without devoting significant computational resources 
to following the evolution of a tenuous corona at large \z\. 
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Equating the right-hand sides of equations (|12p and (|15|) and inserting the expression from equation 
(|13p . this implies /?£ = cr(l — e ±i7rf ) for the isolated sheet, so that the density in real space is obtained 



by taking an inverse transform: 



1 2N ^' 1 

p{z) = olvT ^ 2mt ~^P^ (16) 



2N Z 

z £=0 



■ 2N Z -1 

L £ c-M&oil (17) 



£=0 



2iV z -l 2N Z -1 

^ e z -^:^ e "> (18) 

£=0 £=0 

a<5(z) -a5(z^L z ). (19) 



We see that the first term corresponds to the original density. However, a second term has appeared 
as an image density with the opposite sign from the real (physical) density, located a domain 
length away. This means that to obtain the correct solution for $ on the original domain z = 
(— L z /2, L z /2), we need to prepare twice as large a box in the vertical direction, and implement 
the required image density within the augmented domain, at a distance ±L Z from the physical 
slab. Thus, a density slab at < z < L z /2 would require an image slab at Zimage = z — L z m 
(— L z , —L z /2), and a density slab at —L z /2 < z < would require an image slab at -^image — z ~~~ I j z 
in (L z /2, L z ). Using a similar procedure, we have extended this idea to the three dimensional case 
with an arbitrary density distribution p(x, y, z). The details are described in the App endix. We note 



that the numerical solution agrees with the solution obtained via Green functions (jMiyama et al 



19871 ). and is much faster to compute because only FFTs (no direct sums) are needed. 



2.4. Cooling Function 



To allow for multiphase interstellar gas components, we must solve a thermal energy equation 
th at allows a wide ran g e of c onditions. We use a cooling function for the diffuse ISM derived 
by iKovama &: Inutsukal (|2002l ) , which includes atomic gas cooling for the warm and cold neutral 
medium (WNM, CNM), as well as cold molecular-phase cooling (H2, CO, and dust cooling). We 
include a constant vo lumetric heating rate to represent p hotoelectric heating by diffuse FUV. This 
yields a standard (cf. IField et al.lll969i ; IWolfire et al.lll995l ) thermal equilibrium curve in which there 
is a maximum density and pressure for the warm phase given by 1.0 cm -3 and 5.5 x 10 3 &b cm _3 K, 
and a minimum density and pressure for the cold phase given by 8.7 cm -3 and 1.75 x 10 3 k-g cm _3 K 
(see Fig. BJ. 

HII regions in the real ISM include photoionization of atoms and dissociation of molecules, and 
radiative cooling of photoionized gas and warm molecular gas. These effects depend on chemical 
fractions, as well as dust evaporation. For this work, we are interested primarily in dynamics 
of the neutral media, rather than the details of photoionized gas - including the complexities of 
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ionization front propagation at sub-parsec scales. The main requirement for capturing the large- 
scale dynamical effects of feedback is thus to incorporate photoheating of gas in star-forming regions. 
The simple but expedient approach we have chosen is to expose gas in targeted regions to enhanced 
heating, while simply applying the same cooling function we use for neutral gas. The enhanced 
heating we apply yields thermal equilibrium for the "photoheated" gas at ~ 10,000K (see below), 
which is consistent with the temperatures that would be attained if we had implemented realistic 
(but much more computationally complex and expensive) photoprocesses. 

Cooling and heating timescales often become much shorter than the hydrodynamic time step 
(i.e. the flow or sound crossing time of a grid zone), especially in HII regions, which have a high 
heating rate. For efficiency, we adopt implicit time integration for the heating and cooling operators. 
In a given zone, the integral from the (j) to the (j+1) time step is formally expressed as 

T ^+ 1 C v dT , . 

Atj, (20) 



Tj - nA(T) 

where C v = &b/(t — 1) is the heat capacity per particle. This is a nonlinear equation with respect 
to Tj + i, with Tj and Atj treated as parameters. For this integral, we adopt Simpson's rule and 
solve using the Newton-Rap hson method. 



2.5. Thermal Conduction 



Thermal conduction determines the thickness of interfaces between phases in the ISM, and 
proper incorporation of conduction is essential in numerical simulations of a multiphase medium 



proper incorporation ot conduction is essen tial m numerical simulations or a mu ltipliasc 
which is subject to thermal instability ( Piontek Sz Ostrikerl2004 ; Kovama &: Inutsukall2004 ; 
20081 ). The characteristic length scale set by conduction is the Field length, 



Kim et al. 



Af = 2iTt 



' KT 
n 2 A(T) ' 



(21) 



([Field! Il965l ; iBegelman fc McKed Il990l ). which corresponds to the critical wavelength of thermal 
insta bility. For rea listic values of the conductivity at T ^ 10 4 K, K = 2.5 x 10 3 \/r erg cm -1 s _1 
K _1 ( Parkerl Il953l ). the Field length of 0.19 pc (at density 1 cm -3 and temperature 1,000 K) is 
much smaller than the size of interstellar clouds (~ 1 - 10 pc) and we would need extremely high 
spatial resolution to resolve it - and a correspondingly high computational cost. Instead, we adopt 
an approach somewhat analogous to the use of artificial viscosity (far exceeding the true physical 
viscosity in magnitude) in resolving shocks on a numerical grid. Namely, we adopt a sufficiently 
large numerical conduction coefficient to resolve the Field length on our chosen grid. We choose K 
such that for any simulation with resolution Ax, the Field number N-p = Af/Ax is equal to 1.7 at 
density and temperature typical of thermally unstable gas (we use n = 1 cm -3 , T = 1,000 K). For 
example, the fiducial model Qll has Ax = 1.5 pc and the artificial conductivity gives Af = 2.57 
pc for thermally unstable gas. 
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For low density gas, given our typical values of K the thermal conduction term can become 
greater than cooling/heating terms. In order to limit the conduction in these regions, we adopt 

K 

K = (22) 

1 + n C nt/n 

with n cr ; t = 0.05 cm~ 3 . 



2.6. Stellar Feedback Activity 

The primary focus of this work is to explore the dynamical effects of strong, localized heating 
by OB stars in dense regions of the ISM. Since this heating produces T ~ 10 4 K gas that is 
initially overpressured by a factor ~ 100 or more compared to its surroundings, HII regions expand 
rapidly. This process is believed to be an important source of turbulence both within self-gravitating 
GMCs and in the surrounding diffuse ISM. To study this process, ideally one would implement 
(a) formation of OB stars from dense gas, distributed throughout the space-time domain of the 
simulation; (b) radiative transfer of ionizing photons from the OB stars through the surrounding 
gas, with potentially multiple ionization sources throughout the domain; (c) detailed ionization and 
heating of the gas within HII regions. 

In this first exploration, rather than attempting to model all of these processes in an exact 
fashion, we instead adopt an idealized approach, with the goal of gaining physical understanding. 
First, we apply certain criteria to determine when and where "star zones" on the grid will be turned 
on. Then, we apply strong heating to the gas in the vicinity of each "star zone" for the duration 
of its lifetime. All "star zones" have the same lifetime, i ms , whi ch is set to 3.7 x 10 6 yr, t he typical 



lifetime of OB associations in clouds whose mass is 10 5 Mq ( McKee &: Williams! 1 19971 ). Within 
HII regions, we assume a constant gas heating rate, set via a control parameter Ghii- Each "star 
zone" is therefore essentially a control flag for whether or not strong photoheating is locally applied 
near that zone (which does not move). Rather than solving a radiative transfer problem, we use a 
simple criterion based on the gravitational potential to determine whether gas is subject to strong 
heating. Because our goal is to identify gas localized around star-forming regions, it is necessary to 
subtract out the background gravitational potential and retain just the potential component due 
to an individual self-gravitating cloud. 

The background potential is the potential averaged over horizontal planes. In terms of Fourier 
components, the relative gravitational potential = <1> — ^background is defined as 

<J>«(x) = * fc /kexp(-ik-x). (23) 
Here, /k is the Fourier component of a smoothing window function, 

/(*) = A rr n 1 .. r> (24) 



47 "ol + exp ^(y^T^-ro; 
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where tq is a smoothing length. This window function smooths the HII region within a radius ~ tq. 
Convolution of the relative potential with a smoothing function (or, equivalently, multiplication in 
Fourier space as above) is desirable so that any heating that is applied is resolved on the grid. We 
have adopted tq = 3Ax as providing adequate resolution. 

HII photoheated regions are identified as regions where the relative potential, falls below 
some specified level: 3>W < $hii- We also employ the relative potential for setting one of the 
criteria for turning on feedback: $>W < <3?sf and p > psf must both be met in a given zone for a 
"star zone" to be created at that location. Thus, our recipe ensures that feedback will only occur 
in dense and self-gravitating regions, consistent with the fact that OB stars are observed to form 
only under these conditions. 

For the feedback prescription we have adopted, there are five control parameters: psf, ^SF; 
tms, &mi, Ghii- The detailed estimation of those parameters is described in the remainder of this 
section. 



2.6.1. OB Star Formation Criterion 

We choose a density threshold for star formation as 

p SF = 10 3 /icm- 3 . (25) 

This density is comparable to that of clumps of gas within GMCs. The free fall time at this density, 
1.4 Myr, is sufficiently small compared to the orbital time that structures satisfying this threshold 
evolve independently of the global environment. Note that the local Jeans length 

Aj= f_^Ly /2 = 4.1pcf T ) (26) 

\GpsF J VO.Qkms- 1 / v ; 

must be resolved b y a few zones in orde r to prevent fragmentation occurring as a consequence of 



numerical artifacts (jTruelove et al.lll997l ). 

To obtain an estimate for the potential threshold for star formation, we consider a cloud with 
uniform number density n and radius R c \. For a spherical cloud, the radius is related to the cloud 
mass using 

R* = (^) 1/3 = 19pcx^/X/ 5 3 , (27) 

where the fiducial value of fi2 = n/(10 2 cm -3 ) is chosen using a typical mean density within GMCs, 
and M cl , 5 = M cl /(10 5 M ). 

Since our grid is two-dimensional, the control parameter must be based on a cylindrical 
regions of a given density. For a cylinder of radius R c \, the potential difference between the center 
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and a distance r (> R c \) is ^(O) = 3>(r) — 2TrGpR^{ln(r / R c \) + 1/2]. The logarithmic term cor- 
responds to the potential difference between r and R c \, and the 1/2 is the contribution between 
the cloud's surface and its center. If the cloud is created out of all of the mass originally within a 
disk of surface density £ within a range 2r = L eqv of the cloud center, then nR^p = SL eqv , and a 
fiducial distance for defining the effective zero of the potential is 



irpR : 



Leqv — f] - ^"' 

= 393 P cx^( I;s ^- F ,)" 1 . (29) 

Here, the radius is expressed in terms of that of an equivalent spherical cloud with a given mass. 
If we set the potential at r = L eqv /2 to zero, then the potential at the center of the cloud will be 



$(0) = -$ S F = -2-KGpR, 



In ( ^l] + I 



\2RcJ 



(30) 



For n = 100 cm" 3 and R c \ = 19 pc, the value of $(0) is 3.4 x 10 11 cm 2 s -2 times an order-unity factor 
that varies only logarithmically in the ratio of cloud surface density to mean ISM surface density. We 
choose to adopt a potential threshold for star formation $sf = 3.4 x 10 11 cm 2 s -2 = (5.8km s -1 ) 2 ; 
we have tested sensitivity to the value of $sf, and found that results are insensitive to the exact 
choice, except as described below. 



2.6.2. Definition of Photoheated Regions 



First, consider an HII region centered on the origin of a uniform cloud with spherical crossec- 
tion. If we assume the radius of the HII region is -Rhii (< Rci)> then if the center of the cloud has 
potential ^(0), the potential at the ionization front for a cylindrical cloud with uniform density is 
^hii = ^(0) + 7rG/xR 2 III . (For a spherical cloud, the second term is smaller by a factor 2/3.) Taking 
the difference with $(L eqv /2) in order to represent a relative potential, using equation (|30p . and 
substituting c &(0) — > — $sf (since the criterion for star formation must be satisfied if feedback has 
turned on) this implies that the relative potential at the location of the HII region would be 



$HII = -$SF 



1 



(31) 



Run is given, for example, by the Stromgren radius in a uniform medium: 

3S \ ^ 



-Rhii 



("47rn 2 a( 2 ) 



2.97 pc n 2 27,3 5*49 3 , 



(32) 



where S is the number ionizing photons per unit time and S49 = <S7(10 49 s 1 ), and = 3.09 x 
10~ 13 cm 3 s _1 is the case B hydrogen recombination coefficient at T = 8,000 K ( Spitzerl Il978l ), 
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5*49 i s equal to unity for the ionizing luminosity in a typical 10 5 Mq GMC ( McKee Williams 
19971 ), From equation (|3ip . when the density is comparable to the mean density within a GMC, 
Rmi/Rd < 1, and the HII region is buried deep within the GMC at ||$hh| - ^sfI/'&sf < 1- 

HII regions are initially highly overpressured, however, and will expand rapidly until breaking 
out of the surrounding GMC, creating a blister HII region. For the purposes of considering the 
momentum input to the system, the i?Hll — ► Rci limit is most appropriate for defining the pho- 
toheated region. Thus, we suppose that the HII region has expanded, leaving a very low density 
interior and a shell of radius i?Hii < Rci m which most of the mass has piled up. The potential in 
the interior of the shell is then given by 



HII 



SF 



R 2 
■"•HII 



-"HII 



ln( 



R, 



In 



+ 



(33) 



When -Rhii — * Rci, this expression is of course the same as if we had taken i?Hii/-R c i — ► 1 in equation 
(|3ip . that is, the potential near the surface of the initially-uniform cloud. For convenience, we 
introduce a dimensionless parameter e: 



$hii = -$ S f(1 - e) 



(34) 



where e = (1/2) [In ( -grp I + | ] when iinn — ► Rc\- For the fiducial parameter values discussed 



2R C 

above, e is equal to 0.18. We therefore adopt e = 0.2 as our "standard" value, although we have 
tested how the results differ for much smaller values. 



2. 6. 3. Heating Rate in HII Regions 



During the period that "star zones" are alive, UV photons enhance the heating rate within 
HII regions, defined as described above. For the heating rate in any zone due to UV photons, we 
adopt the on-the-spot approximation given by: 



r 

Gfuv 
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Gfuv erg cm s" 

($« < $ffll), 
(otherwise) , 



(35) 
(36) 



where Gq = 1.0 is the background FUV field in the Galaxy. 



We choose Ghii = 10 3 throughout this paper, although we have also tested results with lower 
and higher Ghii- In practice, the exact value of Ghii is not important, since the purpose of this 
added heating is simply to increase th e maximum density at which a warm phase is present. From 
our Fig. [3] (see also e.g. Fig. 7 of IWolfire et al.1 Il995l ). a value of Ghii = 10 3 boosts this to 
n ~ 10 3 cm -3 . 
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3. Model Parameters 

The models in this paper are characterized by three main parameters: the total gas surface 
density E, the orbital angular velocity f2 in the center of the domain, and the stellar density at 
the midplane p*. The first parameter defines the amount of raw material available for distribution 
between the dense and diffuse ISM phases, while the second two parameters define the galactic 
environment in which the gas evolves in response to the galactic tidal, rotational, and shear stresses. 
The effectiveness of self-gravity in forming massive clouds depends on all of these parameters, as 
well as on the turbulent state that develops as a consequence of star formation feedback. 

The simulation domain we model is a radial-vertical slice through a disk. For the vertical 
direction, we require a domain L z that encloses most of the mass of the ISM, which is confined by 
a combination of stellar gravity and gas self-gravity. The largest scale height (excluding the hot 
ISM, absent in these models) is that of the WNM, and an upper limit is obtained by neglecting 
self-gravity, which yields a Gaussian distribution with scale height: 

/ \ -i/ 2 

Hw = = 95 PC (rknTFr) ( o.i m & pc -3 > (37) 



Here, c w = (kT w / fi) 1 ' 2 is the thermal speed of the warm medium, which typically has T w — 
8, 000 — 10, 000 K. We require L z > 2H W for our numerical models. 

If the ISM consisted only of WNM, then with n w o = Pw,o/p the central volume density, the 
total surface density would be 



£„ = VW*. = M S PC- (^) (^) ^^_ ? )- 1/2 . (38, 



The maximum density for which the warm phase is possible (when Go = 1) is n^ jmax = P W) m&x/ {pc w/ 
1 cm" 3 ; this implies a maximum possible total surface density of warm gas S^ imax = (2G/9*)" 1 / 2 : P Wjmax / 'c w . 
In practice, the midplane density of the warm medium is closer to the value 0.23 cm -3 , which rep- 
resents the warm medium density that is in pressure equilibrium with cold gas at -P c ,min/A;B = 
1.75 x 10 3 cm _3 K. We are interested in multiphase disks similar to those observed in the So- 
lar neighborhood and at smaller radii in normal spirals; hence we choose surface density of at 
least several M Q pc -2 such that the pressure is high enough to permit a cold phase to form, i.e. 
E > (2Gp,)- 1 / 2 P cMn /c w . 

The radial domain should be sufficient to capture the largest-scale gravitati onal instabilit ies, 



which in a disk galaxy are limited by angular momentum. The Toomre length (lToomrelll964l ) is 
the maximum scale for axisymmetric modes in a thin, cold disk; 

4vr 2 GE 4^ 2 GE 

AT = "^ 2_= (4-2g)0 2 ' (39) 
= ° 6 ^25 ^^-1 kpc -i) (lO Mq /pc 2 ) (40) 
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We require L x > At for our numerical models, typically by a factor 1.3. 

The parameters of our models are summarized in Table [TJ In order to cover the three- 
dimensional parameter space and explore environmental dependences systematically, we consider 
four series of models: Q, K, S, and R. For each Series, we hold two quantities fixed and vary a third 
quantity, as follows: 

• Series Q: k/H and y/p^/T, are constant while £ varies; 

• Series K: k and y^/E are constant while £ varies; 

• Series R: k/£ and are constant while £ varies; 

• Series S: £ and are constant while k (and f2) varies. 

The value of the Toomre Q parameter for the gas component, for a radial velocity dispersion 
an, is defined as 

« - («) 

- 1* L "J (^) (^J—,) " ■ (42) 



25km s" 1 kpc" 1 / V7km s" 1 / V1OM pc 

Since Toomre's parameter is proportional to fc/£, Series Q and R would have constant thermal 
Q = kc s /(7tG£) for the gas if the sound speed c s were constant (which is true for warm gas). 
The Q and R series correspond to values of Q = 2.1(ctr/7 km s _1 ). Assuming a constant stellar 
velocity dispersion, £* oc yfpZ, so the stellar Toomre parameter (Q* oc «;/£*) would also have 
the same value for all members of Series Q. As a fiducial model, we choose £ = 10.6 M /pc 2 , 
= 31.2 km s _1 kpc" 1 , and = 0.07 M /pc 3 , similar to conditions in the Solar neighborhood; 
this is denoted as the Qll model in Tabled! 

Relative to conditions similar to those in the Solar neighborhood, we may think of the members 
of Series Q as representing conditions ranging from slightly larger radii down to radii of a few kpc, 
for a disk that has constant Q and Q* - i.e. larger gas and stellar densities at small radii. We 
may think of Series R as models spanning a similar range of radii, except for a disk that has stellar 
density (and the corresponding vertical gravity) independent of radius, while the gas surface density 
increases inward (such that the gas self-gravity can become dominant). We may think of Series 
S as relocation of the Solar neighborhood conditions of gas and stellar density to either further 
in or further out in the galaxy's potential well, where rotation and shear are stronger or weaker, 
respectively. We may think of Series K as choosing a fixed location in the galactic potential well 
(dominated by dark matter), and then varying the gas and stellar surface densities in tandem. 

To initialize our models, we set the density to a uniform value (given by n in Table [1]) through- 
out the domain, and set the pressure to P/k^ = 4, 860 cm _3 K which is in the thermally unstable 
regime. The value of the initial pressure in fact does not matter, since the gas rapidly separates 
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into WNM and CNM due to thermal instability. We also impose on the initial conditions isobaric 
density perturbations (with a flat spectrum at wavenumbers smaller than kL z /2n = 32, and 10 
% total amplitude). The results are also not sensitive to the amplitude or power spectrum of the 
initial perturbation spectrum; this is simply a convenience to seed the initial evolution into thermal 
instability and then Jeans fragmentation of cold gas. In order to reach a quasi-steady state with 
repeated feedback cycles, we run our models for two orbital periods in Series Q, R and K and up 
to tfi na i = 5.57 x 10 s year for Series S. 

3.1. Hydrostatic Models 

Because an important focus of this work is to assess the effects of turbulence, it is important to 
ascertain how our dynamical models differ from the situation in which there are no motions other 
than background sheared rotation. For these comparison models, we calculate the one-dimensional 
hydrostatic equilibria in the vertical direction. These models include heating and cooling, but no 
feedback from star formation. We consider two series, HSP and HSC which have stellar volume 
density either proportional (P) to the square of the gas surface density S or constant (C), 
respectively. Note that Series HSP corresponds to the dynamical Series Q and K, while Series HSC 
corresponds to the dynamical Series R. Details of these model parameters are listed in Table 
Uniform 1024 grids are used for all hydrostatic models. Note that the hydrostatic models only 
allow one-dimensional vertical structure, and require higher spatial resolution than the dynamic 
models because the scale height in the cold layer near the midplane becomes very small when there 
is no feedback. 

4. Evolution of the Fiducial Model 

This section describes the overall evolution of the fiducial model. After the initial thermal in- 
stability, the cold gas collapses into the midplane, due to the vertical gravitational field. The cold, 
high-density midplane gas then fragments rapidly, due to self-gravity. The time scale for this frag- 
mentation is characterized by the Jeans time for a thin disk, tj = c s /GT< ~ 6 Myr(c s /0.3 km s" 1 ). 
The dense fragments collapse, with some of them coagulating, until the feedback criteria are met, 
and heating is turned on to create local HII regions. These HII regions expand, causing the dense 
gas to be swept outward, forming shells that then break up into filaments. New overdense regions 
continue to form, collapse, and be dispersed by feedback. 

Figure Q] shows a snapshot from the fiducial model (Model Qll) at a point after the system 
has reached a quasi-steady state, in terms of the statistical distributions of density, temperature, 
and velocity. The two panels show the temperature and density throughout the domain. The 
contours in the lower panel denote relative potential <1>M: solid and dashed lines show positive and 
negative values, respectively. At the time of this snapshot, there are three large clouds consisting 
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of collections of dense filaments that create minima in the gravitational potential (dashed lines 
in the lower panel). Most of the dense filaments and clumps in the lower panel correspond to 
cold gas in the upper panel. In the upper panel, the orange circle associated with the cloud near 
x = 250 pc shows an active HII region, in which the gas is both warm and dense and hence is 
expanding rapidly. Expansion of shells slows at later times (after the pressure drops to ambient 
values and the enhanced heating turns off). Since most of the time for any given shell is spent near 
the maximum expansion, the widely-expanded structure of the middle cloud is typical, in terms of 
the time-averaged state. 

Figure [2] shows a snapshot of the density and temperature in the same model Qll at a time 
38 Myr later. Overall, the structure is qualitatively similar, although details change because the 
state is highly dynamic. There are still three main collections of filaments; the middle cloud has a 
large shell while the left- and right-side clouds have contracted onto the midplane and have nearly 
reached the point at which new HII regions will be born. 

Three large "clouds" within the 1.16 kpc horizontal length of the domain corresponds to mean 
separations of 390 pc. One might expect the number of cloud entities to be related to the properties 
of star formation feedback, and for our adopted prescription to the parameter e, which effectively 
determines the maximum volume of an HII region: large e corresponds to large HII regions, whereas 
small e corresponds to HII regions only in the immediate vicinity of a potential minimum defined 
by a high-density clump of gas. In a situation with multiple local minima in the gravitational 
potential, if e is large then a single HII region could engulf what would be multiple HII regions in 
the case of small e. Expansion and shell collision of many small HII regions would produce more 
(but smaller) clouds than expansion and collision of a few large HII regions. In fact, when we run 
the same model but set e = 0.02, we find that there are typically 4-5 clouds in the same domain. 
We conclude that the number of large clouds is not very sensitive to e, but since this control 
parameter can only approximately model the effects of real HII regions, the current study cannot 
provide an exact prediction for the size or mass of GMCs. We note that the typical separation is, 
however, in the same range as the two-dimensional Jeans length at the typical velocity dispersion, 
Aj,2D = cr 2 /(GT.) = 340pc(o-/4 km s" 1 ) 2 for this model. 

Figure [3] shows the phase diagram (p — P plane) for the same snapshot in Figure Q] and [2j 
The gray scale is proportional to the fraction of the total mass in the domain that is found at a 
given point in the phase plane. We overlay the thermal equilibrium curves for both the cases of 
"normal" heating (solid curve) and the enhanced heating in HII regions (dashed curve). Clearly, 
most of the gas resides near thermal equilibrium, due to the short cooling times compared to the 
longer hydro dynamical times. 

The range of properties of the gas can also be seen in the Figure HI which shows the probability 
distribution functions (PDFs) of gas density. Solid and dashed lines show mass and volume weighted 
probabilities, respectively. The volume PDF shows that the volume is mainly occupied by warm 
and diffuse gas (WNM) at densities of a few xO.l cm -3 . In terms of the mass PDF, there are two 
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peaks: one corresponds to the WNM, and the other to a cold component at density above 10 cm . 

Figure [5] shows the time evolution of thermal, kinetic and potential energies averaged over the 
domain, for Model Qll. For the potential energy, the background disk potential is subtracted out; 



which correspond to times when new HII regions are born and then rapidly expand. The number 
of spikes corresponds to the number of stellar generations in the model; note that this number must 
be proportional to the domain size. In the second rotational period (i.e., 1 < t < 2 rotation), there 
are 18 generations per 1.16 kpc, or 6 generations per massive cloud per rotation period (i.e. an 
interval of 3.3 x 10 7 yr) if the mean number of clouds is 3 in this fiducial model. 

Other models show similar overall behavior in terms of the evolution of physical structure 
(consisting of clumps and filaments that disperse and re-collect), as well as the distribution of mass 
in the phase plane. As environmental parameters vary along a sequence, however, there are some 
characteristic changes in structure. The most notable difference is that for high gaseous surface 
density cases, clouds are often more physically concentrated (i.e., more compact and dense) because 
of the higher stellar and gaseous gravity. For example, Figure [6] shows a snapshot from Model Q42, 
which has £ and four and 16 times larger, respectively, than the values in the fiducial model 
Qll. The three clouds that are seen in the figure are more compact than in the lower-/)* case. For 
a given $sf and e, increasing the mass of a cloud implies that HII region cannot break out as easily. 



All of our models show a turbulent, multiphase ISM with several generations of feedback from 
photoheating. In this section, we analyze how the statistical properties of those models depend 
on the environmental parameters, both along a given series and from one series to another. The 
statistical properties that we study are based on averages of the fluid variables over space and time. 
First we describe how these averages are defined in general, and then we turn to the particular 
statistics. 



For a variable of A that is averaged over both space and time, we use mass-weighted averages 
defined by: 




There are many sharp spikes in both thermal and kinetic energies 



5. Parameter Dependence of Statistical Properties 



5.1. Space and Time Averaging Procedure 



(A(a)) 



f A dm 

J a 



(43) 





(44) 
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where the argument 'a' denotes a given phase or component of the gas (such as WNM). The time 
averaging is then defined via: 



where the step function gives 1 if M{a) > and 9 = otherwise, so that only intervals in which 
the component is present are included in the averaging. Note that S = t<i — t\ if the component 
'a' is present somewhere in the domain at every interval during the simulation. Thus, we denote 
spatial and temporal averages using angle brackets and overlines, respectively. To avoid the initial 
transients at the beginning of the simulation, we adopt the interval between t\ = ifi na i/2 and 
*2 = ^finai for our time average. For the purpose of averaging, our sampling rate is At = 0.1i ms , 
where t ms = 3.7 x 10 6 yr is the adopted lifetime of "star zone" flags that control photoheating 
feedback. 



The models of this paper focus on dynamics rather than chemistry, so rather than dividing 
the gas into distinct phases, we simply bin it according to density. The neutral gas (i.e. the 
gas that is not within the limits defining HII regions) is separated into four bins. The first bin 
(n < 1 cm" 3 ) corresponds approximately to the WNM, with densities below the maximum for 
which a warm phase is possible in thermal equilibrium. The second bin (1 cm~ 3 < n < 100 cm" 3 ) 
extends up to the maximum density that is in pressure equilibrium with WNM gas, and corresponds 
approximately to the CNM (phase diagrams show that the thermally-unstable regime is not highly 
populated for our models; see e.g. Figs [3] and 0} . The dense medium (hereafter DM) is all the 
gas at n > 100 cm -3 , which in thermal equilibrium is above the maximum pressure for the warm 
phase and therefore only exists in regions that are internally stratified due to gravity. The DM 
gas corresponds approximately to the molecular component of the ISM; we divide it into two bins, 
DM2 (100 < n < 10 3 ) and DM3 (10 3 < n). Gas that is within the limits defined for enhanced 
heating is labeled as ionized gas (hereafter HII). So that the mass fractions f(a) of all components 
add to unity, ^2 a f(a) = 1, we use a slightly different definition from that of equation (145 j) . The 
mass fraction of each component a in (WNM, CNM, DM2, DM3 or HII) is defined as 




(45) 



(46) 



5.2. 



Mass Fractions 




(47) 



t 2 -h 



Figure [7] shows the mass fraction of the various components either as a function of surface 
density (Series Q, R and K) or angular velocity (Series S). For Series Q and R (which most closely 
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correspond to the radial variations found within normal spiral galaxies), at low X the diffuse 
(WNM+CNM) components dominate, while at high S the dense (gravitationally-confined) com- 
ponents (DM2+DM3) dominate. The behavior is somewhat different in Series K, which is highly 
gravitationally unstable and thus extremely active when S is large (since k is constant), leading 
to larger HII and CNM mass fractions at high X. At low X, the behavior in Series K is similar 
to that in Series Q and R. For all series, the HII mass fraction increases at higher X or lower f2, 
corresponding to lower Toomre Q (see Figure [TT1) and hence higher rates of stellar feedback activity. 
The mass fraction of the WNM component secularly declines with increasing X in model series Q, 
K, and R. Even though the models of Series S are most gravitationally unstable at low Q, they 
remain dominated by diffuse gas (CNM) rather than dense gas, because the total surface density 
is relatively modest for this series (X = 15 M pc~ 2 ). 



5.3. Surface Density 

The simulation domain for our two-dimensional models represents a radial-vertical {x — z) 
slice through a galactic disk, such that if we viewed the corresponding galaxy face-on, the surface 
density as a function of radius would be given by S(x) = J dzp(x,z). The area- weighted mean 
surface density in any model is equal to (£)a = J dxdz p/ L x ; this is a conserved quantity for any 
simulation, and is one of the basic model parameters (see column 2 of Table [TJ in general we omit 
the angle brackets and "A" subscript). The value of the surface density weighted by mass rather 
than by area better represents the "typical" surface density of clouds found in the disk. This is 
defined as 

_ fX(x) 2 rfx 

Scloud " /£(*)<**■ (48) 

Figure [8] shows X c i ou d as a function of X for all model series. Interestingly, we find that X c i ouc i does 
not strongly depend on parameters (either X or Q) throughout the four model series. The largest 
value of X c i ou d is 3OOM pc -2 and the smallest is 5OM pc -2 , although for most cases the range is 
even smaller; X c j ouc [ = 70- 15OM pc~ 2 . This factor-of-two range for X c i OU( i is significantly smaller 
than the factor-of-six range of mean surface densities, X = 7.5 — 42M pc -2 . The range of X c i ou d 
is also similar to the typical observed surface densities of giant molecular clouds (see discussion in 

ED- 

This weak variation of X c i ouc j among the various series suggests that it is star formation feed- 
back, rather than the feedback-independent parameters, that determines the typical surface density 
of clouds. In particular, other tests we have performed suggest that it is the gravitational potential 
threshold for star formation $gp that most influences X c i ou d. A model in which <I>sf — ► (with 
any e) would have photoheating events independent of the local ISM properties; the consequent 
expansion of HII regions would thoroughly mix gas so that X c i ou d — > X. This is indeed what we 
find when we run models with |3>sf| a factor 10 below our adopted value. On the other hand, larger 
values of |3>sf| require more massive and compact clouds in order to have star formation, which 
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would raise E c i ud- Tests with |<3?sf| a factor 10 above our adopted value indeed result in larger 
^cioud (although only by a factor ~ 2). The dependence on e is much weaker than the dependence 
on $sf; reducing e by a factor 10 changes £ c loud by only tens of percent, at our standard ^sf- 

The comparison between Series Q and Series R is also interesting, in this regard. The difference 
between these two series is that the stellar density increases with E in Series Q, while p* 
is constant throughout Series R. Based on the larger £ c i ou d value for the largest £ in Series Q 
compared to Series R, when the stellar density is increased, the surface density required in order 
to form clouds also increases. 

5.4. Temperatures 

Figure [9] shows the space-time averages of temperature for the components we have defined via 
density bins, (T(a)), where the argument 'a' denotes WNM, CNM, DM2, DM3 and HII. Throughout 
the model series, the temperatures for the most dense and most diffuse components are fairly 
constant; we find T = 6, 000 - 7, 000 K for WNM, T = 20 - 40 K for DM2, and T = 10 - 20 K for 
DM3. For the CNM component (which in fact includes thermally- unstable gas when it exists), the 
range is somewhat larger, T = 100 — 400 K, reflecting the larger range of conditions for this gas. 
The gas which is subject to enhanced heating has mean temperatures for most models of 4,000 - 
8,000 K. 

The link between density and temperature in our models implies that the components we have 
defined via density bins also approximately correspond to natural ISM phases. This is because 
much of the gas mass is close to thermal equilibrium (see Figure [3| , and we have chosen the bin 
edges so as to match up to points in the phase plane with physical significance. Temperature PDFs 
that we have constructed show a bimodal distribution, as is expected based on the cooling function. 

5.5. Turbulent Velocities 

Expansion of HII regions feeds kinetic energy into the ISM. This kinetic energy is not imparted 
solely to expanding HII bubbles and shells surrounding them, but is shared throughout the ISM 
as turbulence. Our models provide a first look at the results of this form of turbulent driving. It 
is interesting to examine how the turbulent amplitudes vary from one component to another in 
a given model, and how the overall levels vary between models with different feedback rates as a 
consequence of different system parameters. 

Figure [10] shows the turbulent velocity dispersions for all series, defined for each component 

as: 

v(a) = ^(v x {aY + v y (aY + v z {a) 2 ), (49) 
where the argument 'a' denotes WNM, CNM, DM2, DM3 and HII, and v y = v y + fix in order 



-22- 



to subtract out the velocity of unperturbed (sheared) rotation about the galactic center. The 
azimuthal velocities are excited by Coriolis forces so that the relation for epicyclic motions 

v y k 1 

V x - 20 = 7! (50) 



should apply (jBinnev <fe Tremaind 119871 ). and we have checked that this is in fact satisfied. We 
note that the velocity dispersion for each component is computed by summing over the whole 
domain. Thus, the measured velocity dispersions are larger than they would be within smaller- 
scale clouds in the system. However, we have found that there are not large contributions to the 
velocity dispersion from velocity differences of widely-separated regions; this is because turbulence 
is driven by the expanding HII regions, such that the maximum correlation scale is comparable to 
the effective thickness of the disk. For example, if we divide the domain in Fig. Q] horizontally into 
eight equal parts, the mean velocity dispersion of all gas within these sub-domains is ~ 98% of that 
of the domain as a whole. Considering just the dense gas, the velocity dispersion for subdomains 
is ~ 75% of that in dense gas for the whole domain. For the three large clouds seen in Fig. [6l the 
mean internal velocity dispersions are an order of magnitude larger than the dispersion in mean 
velocities. 

In general, the densest component (DM3) has the lowest velocity dispersion, with the next- 
densest (DM2) the next-lowest. The value of the velocity dispersions for the dense components are 
highly supersonic, and are similar to (or slightly below) those that are observed within real GMCs 
(see JO). The CNM component in our models typically has higher turbulence levels than the WNM 
component, because the former is in closer (space-time) contact with energy sources. Because 
turbulent motions in our models are driven by the pressure of photoheated gas, P = p<? s ~ pv 2 , the 
turbulent velocities have an upper limit of the sound speed in gas heated to 8,000 K, chii ~ 7 km s . 
Since the driving is intermittent, this upper limit is not usually reached; mean values for the diffuse 
gas are closer to ~ 5 km s . The diffuse-gas velocity dispersions in our models are lower by about 
50% compared to observed levels, indicating (consistent with expectations) that other turbulence 
sources are important in the real diffuse ISM. 

The model series Q and R, which have E oc f! (and thus effectively constant gaseous Toomre 
Q if the velocity dispersion is constant) show velocity dispersions that are insensitive to the value 
of X. Series K, on the other hand, has much higher turbulence levels for large S. This is because, 
with constant Q (= k/\/2) the high-S models are quite susceptible to gravitational instability (in 
terms of Q); this leads to very active feedback which then raises the velocity dispersion. A similar 
physical effect is seen in Series S: the velocity dispersion is highest at low f2, since these are the 
most gravitationally-susceptible models among the series. We discuss measurements of the effective 
Toomre Q values that account for turbulence, in the next subsection. 
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5.6. Effective Toomre Q Parameters 



For a rotating disk that contains only thermal pressure, susceptibility to growth of self- 
gravitating perturbations depends on the Toomre Q parameter, defined by setting an equal to 
the thermal sound speed in equation (|42p . An infinitesimally-thin gas di sk is unstable to ax- 
isymmetric perturbations if the value of this thermal Q-parameter is < 1 (jToomrd 1 19641 ) . Non- 



itv (Goldreich &; Lvnden-Bell 


1965b 




Joe & Solomon 


1984: Rafikov 2001: Kim & Ostriker 2001: 


Kim et al. 


2002, 


2003; 


Kim & Ostriker 


2007; 


Li et al. 


2005) allowing growth at higher 0, while 


nonzero disk thickness suppresses gravitational instability ( 


Goldreich & Lvnden-Bell 1965a: Kim et al. 



2002; iKim &: Ostriker! 120071 ). lowering the the critical Q value. Allowing for all of these effects, 
threshold levels measured from simulations are Q ~ 1.5. 

Turbulence at scales below the wavelength of gravitational instability can also help to suppress 
the growth of large-scale density perturbations, by contributing to the effective pressure. Since 
the original Toomre parameter is arrived at based on effects of radial pressure gradients, only 
the radial component of the velocity dispersion should be added to the thermal velocity disper- 
sion in defining an effective Q (see eq. 02]). It is natural to expect galactic disks to self-regulate 
the values of the effective Q: growth of self-gravitating instabilities subsequently leads to star 
formation and energetic stellar feedback, which drives turbulence, raises Q, and tends to sup- 
press further GMC formation. Indeed, the suggestion that galactic star formation is self-regulated 
through t urbulent feedback dates back to the ea rliest work o n large-scale instabilities in galactic 
gas disks (jGoldreich Lvnden-Belllll965bl ). with iQuirkl (|1972l ) making the related suggestion that 
galaxies deplete their gas until they reach marginal stability. The self-regula tion processes are 
complex, but they hav e begun to be studied in recent numerical simulations (e.g lWada et al.ll2002l : 
Tasker fc Brvanlhoosl ). 



We use the results of our models to measure the values of the effective Toomre parameter in 
the saturated state. We compare four different measurements of Q in each model. The first is 
closest to Toomre's original definition for a gaseous medium in that it is based on thermal velocity 
dispersion; since our medium has components at differing temperatures, we use a mass-weighted 
thermal velocity dispersion: 

K 



Q T (total) 



' 7fcB (T). 



(51) 



7rGs y /i 

The second measurement incorporates turbulence, again including all gas and weighting by mass: 



Q (total) 



/i 



CO + <4>- 



(52) 



For the third and fourth measurements, we consider only the dense gas for both the numerator and 
denominator: 



Q T (n > 100) 



7 fc B 



(T(n > 100)) 



vrGS f(n > 100) 



(53) 
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and 

/ 7^ 



f-(T(n > 100)) + (v 2 R (n > 100)) 
Q(n > m = __-L^ __ . (54) 

Here, the mass fraction of dense gas f(n > 100) is given in eqn (|47p . The turbulent velocities 
dominate the dense gas when they are included, since the thermal sound speed is < 0.5 km s -1 
whereas turbulent velocities are several times larger; see Fig. [TUl Note that k and X are constant 
in time for any simulation. 

Figure [TT1 shows the measured value of these four quantities Qt (total), Q (total), Qt("- > 100), 
and Q(n > 100), for all of our models. In general, we find that the saturated-state values when 
turbulence is included are near unity. The only significantly larger values are for the low-E models 
in Series K, which have large k, and hence the thermal value QT(total) is large; when turbulence is 
included this is raised even more. Since Series Q and R have k/E constant, the value of Qt (total) 
is simply proportional to the mass-weighted thermal velocity dispersion. The increased fraction of 
cold gas at high X leads to a corresponding decrease in Qt (total) . Since the thermal and turbulent 
velocity dispersions of the dense gas are small compared to those of the diffuse components, the 
dense components contribute to Qt (total) and Q (total) mostly by lowering the mass fraction of 
the diffuse components, in the numerator. Because turbulent contributions are positive-definite, 
Q (total) > Qt (total). 

The strongest evidence of self-regulation by feedback-driven turbulence is seen in the saturated- 
state results for Q(n > 100). With the low values of the temperature in the dense component, the 
thermal-only values for the dense gas are mostly Qt(^ > 100) < 0.5. When turbulence is included, 
however, the saturated-state value of Q(n > 100) is between 1 and 2 for almost all models. This is 
consistent with expectations for marginal instability. We note in particular that velocity dispersions 
in Series K (see Fig. [10]) vary strongly with £ (by a factor ~ 5), while Q(n > 100) varies weakly 
with £ (by a factor < 2); feedback evidently self-adjusts in these models so as to maintain a state 
of marginal gravitational instability. 



5.7. Virial Ratios 



In a self-gravitating system that approaches a statistical steady state, the Virial Theorem 
predicts that the specific kinetic and gravitational energies E and W will be related by 2E + W = 0; 



this i s modified when magneti c terms are present (jChandrasekhar &: Fermi 



1956 



1953 



Mestel & Spitzer 



McKee &: Zweibell 1 19921 ). Classically, the Virial Theorem has often been assumed to hold 



within individual GMCs in order to obtain estimates of their masses, and indeed th is yields values 



Solomon et al 



that are consistent (within a factor ~ 2) with other measures of the mass (e.g. 
19871 ). If individual GMCs are short-lived, however, they may not satisfy the Virial Theorem 
because the moment of inertia tensor is changing rapidly enough, and/or surface terms are large 
enough, to be comparable to the kinetic and gravitational energy integrated over the cloud volume 
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(jBallesteros-Paredes et al.lll999l : iMcKee fc Ostrikerll2007l ; iDib et al.ll2007l ). When averaged over an 



ensemble of clouds, there will be (partial) cancellation of surface and time-dependent terms, as 
they appear with opposite signs for forming and dispersing clouds. An added complication is that 
self-gravitating GMCs form out of diffuse gas, and when they are destroyed (whether after a short 
or long time) they return to diffuse gas; thus the different terms in the Virial Theorem may be 
observed in different tracers depending on whether diffuse gas is primarily atomic or molecular. 

Here, we consider virial ratios 

IE 

n ^W\ (55) 

separately for each component of the gas in our models. The term E includes both thermal and 
bulk kinetic energy, computed via a space-time average as: 

1_ _ 1 



E = -(«! + vl + vl) + -(P/p), (56) 

2 y 7 — 1 

and for W only the perturbed gravitational potential is used in computing the space-time averaged 
value of the energy: 

W = 1(503). (57) 

As for the other statistical properties we have considered, we measure 1Z separately for each compo- 
nent (separated into density bins) of the system. Figure [12] shows the virial ratio of each component, 
for all models in all series. Note that 1Z < 2 and 1Z > 2 imply gravitationally bound and unbound 
states, respectively for any component. As we do not separate the contributions to the potential 
from the different density ranges, a given component may be bound within a potential well that 
is created by more than one component. Strictly speaking, the factor 1/2 in equation (|57p applies 
only for self-potentials. 

As expected, the lowest-density WNM component (n < 1 cm -3 ) has 1Z very large (above 
100), and the intermediate-density CNM component (1 cm~ 3 < n < 100 cm -3 ) also is non-self- 
gravitating, with 1Z in the range ~ 5 — 10. The HII (photoheated) component generally has values of 
1Z ~ 10, similar to that of the CNM component. Although it has very large thermal energy (much 
greater than the CNM), the photoheated gas by definition resides within deep parts of the gravi- 
tational potential well. The two dense components, DM2 (10 2 cm -3 < n < 10 3 cm -3 ) and DM3 
(10 3 cm -3 < n), on the other hand, are consistent with being marginally or strongly gravitationally 
bound, with 1Z 2 and 1Z ^ 1, respectively. For the majority of models, the value of 1Z for the 
densest component, DM3, is quite near unity, indicating consistency with virial equilibrium for the 
component as a whole. For a few models, 1Z is as low as 0.3 for the DM3 component; this indicates 
that the dense gas is transient, with dense regions being rapidly dispersed into lower-density gas 
by the feedback process. Overall, we find no significant differences in the trends for 1Z between 
different series or different models within any series. There is a weak correlation between Q and 
1Z, with lower-Q (more unstable) models having slightly higher virial ratios. 
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5.8. Vertically- Averaged Density and Free-Fall Time 

Although the ISM consists of many phases at different densities, all of this gas resides within a 
common potential well which tends to confine material near the galactic midplane. The scale height 
of each phase depends on the support provided by thermal and kinetic press ure (plus support by 



magn etic stresses and cosmic rays, although these may be less significant). In lKoyama &: Ostrikei 



(|2008l ) we consider in detail the vertical distribution of gas within our models, and show that 
vertical equilibrium is a good approximation for the system as a whole, provided that appropriate 
accounting is made for the differing velocity dispersions of different components. We also discuss 
dependence of the mean scale height on model parameters. 

For the purpose of assessing gravitational timescales of the overall ISM system, it is useful to 
measure the density when averaged over large scales (i.e. a volume at least comparable to the scale 
height). To evaluate this volume average in our models, we first compute the vertical scale height, 
defined using the following averaging: 



#ave = Z ° neSpZ2 (58) 

y zones P 

where z is the vertical coordinate relative to the midplane. We further average the values of H mc 
over time. For a Gaussian density profile, p(z) = poexjp(—z 2 /2H 2 ), the midplane density is related 
to the surface density and scale height by po = ^/(V^ttH), and the mass- weighted mean value of 
the average density is given by po/V2. We therefore define an average density in our models as: 

S 

Pave = ,- (59) 



(see also Appendix in lKoyama &: Ostrikeii (120081 )). 



Figure [13] shows the vertically-averaged density for all models in all series. In general, we 
find that the average density increases with the total surface density of gas in the disk. A slightly 
shallower increase of p ave with £ is obtained in Series K compared to Series Q, which can be 
attributed to the large velocity dispersions in strongly unstable (small Q) disk models. Series R also 
has a shallower slope than in Series Q, because the stellar gravity does not increase at large S in the 
former. For reference, we also plot in Figure[13]the values of the vertically-averaged density from our 
comparison hydrostatic model series. The slope of Series HSC (lower-left panel) is shallower than 
that in Series HSP (top panels), again because the stellar density does not increasingly compress the 
gas at large £ in Series HSC. The volume-averaged densities of the dynamic models are lower than 
those of the hydrostatic models by up to an order of magnitude; the differen ce increases at large 



surfac e density where turbulence plays an increasingly important role (see also iKovama fe Ostrikei 



(12008 )) ■ 

Using the mean density and the definition of the free-fall time, 



" (jy V2 - <60> 
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we can calculate the free-fall time for the system as a whole, ifj(p a ve)- Since p ave increases with 
£ in our models, iff(p ave ) will decrease with increasing £. Because star formation requires gas to 
become self-gravitating, a widespread notion is that the star formation timescale, when averaged 
over large scales in a galaxy, will be proportional to the large-scale average of tg, i.e. iff(p ave )- 
Since star formation takes place within GMCs that have much higher density than the mean value 
in the ISM, the conditions that control star formation where it actually takes place are not those of 
the large-scale ISM. Thus, implicit in the notion that star formation times should be related to the 
large-scale mean iff(p ave ) is the idea that the formation of GMCs (on timescales closer to iff(p a ve)) 
is the principal means of regulating star formation. If the star formation efficiency per GMC is 
constant, then the GMC formation rate would control the star formation rate. Alternatively, the 
star formation rate might be related to the large-scale iff(p avc ) if the densities within GMCs are 
proportional to the large-scale mean densities of the ISM, (pgmc) °c Pave- 

Another important dynamical timescale in disk galaxies is the orbital time, t OT ^, = 2tt/0,. 
Growth of large-scale self-gravitating perturbations in disks in fact occurs at timescales longer 
than iff(/O a ve) (provided pressure limits small-scale collapse), and more comparable to t or b = 27r/S7. 



Observations (jKennicuttl 1 1 998bl ) show that empirically-measured star formation timescales in disk 
galaxies tend to be correlated with the orbital time, with ~ 10% of gas being converted to stars per 
galactic orbit. It is useful to compare tff(p ave ) with i or b in our models. Figure [TH shows the ratio of 
%(Pave)Aorb for all hydrodynamic and hydrostatic series. For the hydrodynamic series, the typical 
ratio is 0.06 — 0.2; for the hydrostatic models, the densities are much higher at large £ so that 
%(Pave)Aorb ~ 0.02 — 0.2. For series Q and R, the ratio iff(p av e)Aorb varies relatively weakly with 
£, and lies in the range 0.1 — 0.2. The comparison hydrostatic models for these series also show 
*ff(Pave)Aorb varying only modestly with £. For these series, i or b oc l/£. Since turbulent velocity 
dispersions do not depend strongly on £ for Series Q and R, p ave does not strongly depart from 
a scaling oc £ 2 , yielding behavior similar to tg oc £ oc t or b- Interestingly, the K series, which has 
constant i or b, shows a smaller range of p ave than the Q series. This is indicative of self-regulation: 
high feedback activity in the highest-£ models of series K yield high turbulent amplitudes, which 
lead to lower values of p ave . As a consequence, the ratios of iff/i rb are more modulated in the 
hydrodynamic models for Series K than in the corresponding hydrostatic series. 



6. Implications for Star Formation 

In the present work, we do not explicitly follow star formation. Nevertheless, it is interesting to 
explore the consequences of our statistical results, within the context of recipes that are commonly 
adopted for star formation in numerical models. We compare estimates of the implied star formation 
timescale both to observations and to various fiducial dynamical times. 
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6.1. Star Formation Rates and Timescales 

A common practice in numerical simulations of galactic evolution is to assume that the star 
formation rate per unit volume (in a computational region) is proportional to the gas density per 
unit volume divided by the free-fall time at that density. When a minimum density threshold for 
star formation is imposed, the total star formation rate (SFR, in mass of new stars per unit time) 
takes the form 

M. = eff(fth)M(P > fth) (61) 
iff (Ah) 

provided that the density PDF decreases above the threshold, so that most of the star forming 
activity is in gas near p t ^. Here, the star formation efficiency per free-fall time, £ff(Ah)> is an 
arbitrary constant parameter that is adopted, generally by comparing to observations. In practice, 
the parameter eg (Ah) m this sort of recipe enfolds many different effects that limit star formation 
compared to the fastest possible rate. Within GMCs, turbulence and magnetic fields limit the rate 
of core formation and collapse, and feedback from star formation limits GMC lifetimes; at larger 
scales, dynamical processes in the diffuse ISM limit GMC formation. Depending on the value of 
the threshold density, either more (low pth) or fewer (high pth) processes are implicitly packaged in 
the single efficiency parameter eg(/9th)- 

For a given star formation rate, the star formation timescale is defined by dividing the total 
gas mass by the total SFR, M*: 

M t Qt _ M t Qt tff(pth) 
SF " M, " M(p> Pth )e s (p th y 1 j 

1 *ff(Pth) _ r S F ^ 



f(p > Ah) eff(Ah) eff(fth)' 

The latter expression uses the mass fraction / as defined in equation (|37|) : M tot is the total gas 
mass. Because / and eg- are, by definition, less than 1, the star formation time always exceeds 
the free-fall time at the threshold density. Since the efficiency per free-fall time is arbitrary (from 
the point of view of simulations), it is convenient to introduce rgp = cs(pth)tSF such that rgp = 
iff(pth)//(p > Ah)- This scaled star formation time then depends only on the choice of density 
threshold and the fraction of the total gas mass above this threshold. 

In numerical simulations, the density threshold for star formation is an arbitrary parameter; 
what difference does the choice of this value make to the resulting SFR? To address this question, 
we first compare values of rgp using two different thresholds, n = 10 2 cm 3 and n = 10 3 cm 3 . 
Both threshold values are large enough that the gas at these densities is in gravitationally-bound 
structures, based on the results shown in Fig. [T2J Figure [15] shows the values of the scaled star 
formation time, tsf, for all models. For our chosen density thresholds, tsf is in the range 3x 10 6 — 10 7 
yr for all models. The true star formation time, isF> exceeds tsf by a factor 1 ; the value of eg(/?th) 
must then be quite small (< 0.01) for igF to be > 10 9 yr. Also, since tsf is larger for the threshold 
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choice n = 10 2 cm 3 than n = 10 3 cm 3 , the value of eff(pth) would have to be smaller for the higher 
density threshold, in order to yield the same value of isF at a given E. Note, however, that while 
the thresholds differ by a factor 10, the values of tsf (and hence required eff(pth)) differ by less 
than a factor 2. This reflects the fact that / decreases with increasing n; between n = 10 2 cm 3 
and n = 10 3 cm 3 , our results imply a dependence /(> n) oc n _s with the range of s = 0.2 — 0.5. 
Alternatively, we can think of our results requiring a choice for eff(pth) oc with the range of 
r = — 0.3 in order for the SFR to be independent of the choice of threshold at high densities. 

Other aspects of the results shown in Figure [15] are also interesting. First, it is evident that 
tsf depends only weakly on both surface density (Series Q, K, R) and angular velocity (Series S). 
For Series Q and S, tsf decreases with increasing S. Interestingly, the hydrostatic models show a 
similar range of tsf to the dynamic, turbulent models. The fact that tsf is not strongly sensitive 
to environmental conditions (total available gas content, local shear rate, level of turbulence, etc.) 
may help to explain why empirical SFRs show such a regular character in observed galaxies, in spite 
of widely- varying local conditions. Conversely, the insensitivity of tsf to conditions within a model 
has implications for evaluating theoretical results: successfully reproducing levels of star formation 
similar to observations is not a critical discriminant of how well a simulated galaxy resembles a 
real system. Our hydrostatic models bear minimal resemblance to real galaxies, yet for a choice 
of eff(pth) ~ 0.01 consistent with observed efficiencies in CO-emitting gas in GMCs (which have 
densities in the range n = 10 2 — 10 3 cm 3 ), the resulting star formation times are ~ 4 x 10 8 - 10 9 
yr, similar to the observed range of isF for S comparable to the range in our models. 

To connect more directly to the way observed SFRs are normally presented, in Figure [16] 
we show results for scaled surface density of star formation as a function of surface density of 
gas (Series Q, K, R) and angular velocity (Series S). The scaled SFR per unit area is defined as 
EsFRAff(/0th) = Etot/fSF) where the SFR is taken to follow 

v E to t £ff(Pth)Etot/(p > Pth) 

Esfr = - — = —, r • (64) 

*SF iff (Pth ) 

As before, we compare results based on two different threshold density choices, and also show 
the results from hydrostatic models. Observations are typically fitted to power laws of the form 
Esfr oc E 1+p . For reference, we show slopes with 1 + p = 1 and 1.5. For each model series 
and each value of p t h, we fit a power-law index. We find 1 +p equal to 1.32, 1.43 (Series Q for 
n = 10 2 , 10 3 cm 3 ), 0.94 (Series K for n = 10 2 cm 3 ), and 1.24, 1.19 (Series R for n = 10 2 , 10 3 cm 3 ). 
For the hydrostatic cases, the indices are 1.38 (Series HSP) and 1.39 (Series HSC) at n = 10 2 cm 3 . 
As we shall discuss further in $7] these results are similar to the observed ranges of power-law 
indices that have been reported. We note that Series Q and R show more regular behavior than 
Series K. This reflects the different environmental parameters that are inputs to the models: in 
Series K, the epicyclic frequency is held constant, while in series Q and R we adopt a scaling f2 oc E. 
For eg ~ 0.001 — 0.01, both the magnitude and scaling of the Esfr vs. E results in Series Q and 
R are similar to observations. 
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6.2. Comparison of Timescales 



In we investigated the relationship between the mean large-scale surface density £ and 
the star formation time based on the amount of high-density gas (at n > 10 2 cm -3 within a zone). 
This gas may be considered immediately eligible for star formation, since it is cold and found in self- 
gravitating systems. As noted in §5.81 if formation of massive, cold, gravitationally bound systems 
is the principal throttle for star formation, then star formation times would also be expected to 
vary with the timescales for GMC formation. 

GMC formation is a complex process, and to date no simple formula has been obtained for 
the formation rate. Instead, several different "large-scale" dynamical times are commonly invoked 
to obtain estimates of the GMC formation time. These include the free-fall time at the large-scale 
mean density, ts(p a . ve ), the Jeans time based on the surface density and the gas velocity dispersion, 
ij, and the orbital time i or b, which is generally related to the epicyclic and shear times. It is 
interesting to explore how our measurement of tsf compares to each of these times, as a function 
of the independent parameter in each series. 

We begin with the orbital time, i or b = 2tt/Q. Figure [T7] shows the ratio between tsf = 
€s(pth)tsF (for the two different density thresholds p t ^) and the orbital time. In Series K, E is the 
independent parameter, but t or b is independent of £; thus the ratio is simply a rescaled version 
of tsf shown in Fig. [T5l In Series Q and R, the independent parameter is £, and i or b oc £ , 
so TsF/^orb oc tsfE. In Series S, i or b is the independent parameter. Although the variation with 
the independent parameter is moderate in all the series, the ratio is not constant, and for some 
series shows secular trends. Namely, for Series Q and R, which showed a trend of decreasing tsf at 
increasing E, rsF/^orb increases at larger E. Thus, assuming that the star formation time is oc i rb 
would increasingly overestimate the true star formation rate (presumed to depend on the amount 
of dense gas) as E increases. 

We next consider the Jeans time for a disk, tj = cr/(G£), where a is either the thermal or the 



total (radial) turbulent velocity dispersion, c s or y c 2 + v^. We note that the ratio tj/t rb is given 
by Q/(2y/2), where the Toomre parameter Q is either based on mean sound speed or the total 
velocity dispersion (see §5.6|) . Figure [181 shows the ratio between tsf and the Jeans time, using 
the total velocity dispersion. Again, strong secular trends with E are evident; tj is not a good 
predictor of the star formation time. 

Finally, in Figure [19] we show the ratio between tsf and and the free fall time at the vertically- 
averaged large-scale mean density ( §5.8{ see Fig. [T¥|) . Although the values of this ratio are closer 
to unity than rsF/^orb and tsfAj, we still see that tsf / 'tsiPave) is not constant as a function of E. 
When we compare £/%(p a ve) to the scaled star formation rates based on high density gas shown 
in Fig. (|16p . we find a steeper rise with E, close to oc E 2 in Series Q and R and slightly shallower 
in Series K. 




If star formation is regulated by the collection of diffuse gas into self-gravitating regions, then 
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strictly speaking one would expect specific star formation rates to vary proportional to the fraction 
of diffuse gas divided by the GMC formation time (estimated just including the diffuse gas). The 
above comparisons between tsf and timescale estimates based on mean large-scale properties can 
be corrected to account for this, yielding the ratios t S f(1 - /dense )/*orb, t S f(1 - /dcnsc) 3//2 Aff (/Wc), 
and tsf(1 — /densc) 2 /*j- We find, however, that these ratios also are non-constant in any series, 
although the correction factor does tend to flatten out the secular rise with increasing S in the 
comparison to the free-fall time. 

Taking all of our results together, we conclude that several commonly-used estimates for galac- 
tic star formation timescales based on large-scale mean galactic properties may have only limited 
utility for making detailed predictions of star formation rates. That is, the orbital time, the Jeans 
time, and the free-fall time based on the vertically-averaged density are not proportional to the star 
formation time based on the amount of dense, gravitationally bound gas that is present. Simula- 
tions with insufficient resolution or limited physics may therefore not be able to provide accurate 
predictions of star formation rates, if they do not capture processes at small enough scales to 
represent dense, gravitationally-bound structures. 

7. Summary and Discussion 

We have developed a numerical hydrodynamic code to study the life-cycle of multiphase, 
turbulent interstellar gas in disk galaxies; our model includes gas self-gravity, the vertical gravity 
associated with a fixed stellar disk, radiative cooling and heating in the temperature range of 
10 K < T < 10, 000 K, sheared rotation in the background galactic potential (we adopt a flat 
galactic rotation curve), and a prescription for feedback in the form of HII regions that originate 
within massive, dense clouds. Our simulation domains represent slices in the radial-vertical plane of 
a galactic disk. We focus on scales large enough to include vertical stratification {L z up to 410 pc) 
and significant shear of the disk angular velocity (L x up to 1.6 kpc), but small enough to resolve 
sub-structure within dense, self-gravitating clouds that form (typical zone resolution is ~ lpc). 
Our models are 2.5-dimensional, in the sense that all three components of the velocity are time- 
dependent functions. For feedback to occur, we impose thresholds on both the local volume density 
and on the gravitational potential, so that HII regions only occur within massive clouds (consistent 
with observations). The expansion of HII regions drives turbulence in all the components of the 
gas. We have performed a large suite of simulations, covering a factor of six in gas surface density 
S. In order to explore the dependence of ISM properties on galactic environment (in particular, 
the stellar vertical gravity and the angular momentum content of the gas), we have considered four 
different model series. In Series Q, we vary S, stellar volume density p*, and the disk rotation rate 
f2 in tandem. In Series K, we vary £ and p* together while holding Q fixed. In Series R, we vary £ 
and f2 together while holding p* fixed. Finally, in Series S the values of £ and p* are held constant 
while Q is varied. 

Our main conclusions, and their relation to other recent work, are as follows: 
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1. Density, temperature, and pressure distributions 



We find that in spite of time-dependent effects, the d ensity and t e mper ature distributions of 
the gas retain bimodal profiles reminiscent of the classical Field et al.l (|1969l ) two-phase model of 
the ISM. Although large-amplitude turbulence heats and cools via PdV work and entropy produc- 
tion in shocks, most of the gas (by mass) remains near the curve in the pressure-density phase 
plane that is defined by radiative equilibrium: n 2 A — nT = 0. This is possible because the cooling 
time is generally shorter than the turbulent dynamical times, for our models. If turbulent compres- 
sions or expansions were more extreme in magnitude and also rapid in time compared to radiative 
times, then these adiabatic changes would le ad to density-temperat u re pa ir s that more strongly 



depar te d from thermal equ ilibrium (see e.g. ISanchez-Salcedo et al.l (|2002l ); lAudit k, Hennebelle 



(|2005l ) ; iHeitsch et al.l (120081) f o r a di sc ussion of the depen d ence on var ious physical timesc a les in - 
volved). lAudit fc Hennebelle! ifcoOfJ ). IPiontek fc Ostrikeri ifcoOfJ ). and IPiontek fc Ostrikeri (J2003) 
similarly found that even with large-amplitude turbulence - and no direct heating of the gas 
- mass- weighted density and temperature PDFs have two (broadened) peaks, altho u gh for very 



high a mplitude turbulence the trough tends to fill in (e.g. iHennebelle &i Auditl 120071 ). iKim et al 
(|2008l ) have emphasized the importance of maintaining sufficient numerical resolution, as numer- 
ical diffusion associated with flow across the grid broadens warm/cold interfaces, populating the 
thermally-unstable regime in the phase diagram. We note that bimodal character is most easily 
seen in mass- weighted rather than volume- weighted PDFs, although many of the results in the lit- 
erature show only volume-weighted PDFs. Our results on the bimodal thermal distribution of gas 
are c onsistent with observations of atomic HI in the midplane of the Milky Way (IHeiles fc Troland 
2003), which at the same time show interesting evidence of out-of-equilibrium gas, particularly at 
high latitudes. 

2. Turbulence 

We find that appreciable turbulence can be excited in all components of the gas. The values 
of the velocity dispersion in dense gas [n > 10 2 cm~ 3 ) of ~ 2 — 4 km s _1 are sim ilar to, though 



slightly lower than, those observed in Milky Way and other lo cal-group GMCs (e.g. ISolomon et al 



1987 



Sheth et al 



2008 



Bolatto et al. 



2008 



Heyer et al. 2008). For lower-density gas in our models, 



velocity dispersions are slightly higher, but still lower by a factor two compared to observed velocity 
dispersions of ~ 7 — 10 km s -1 of Solar-neighborhood warm an d cold HI seen in 21 cm emission 
and absorption (e.g.. IHeiles Trolandll2003l ; iMohan et al.ll2004l ). 



It is not surprising that the turbulence levels in our simulations are only moderate, given that 
we have included only one of the many sources o f turbulence that is presen t in the ISM. Turbu- 
lent d riving in the ISM has been reviewed by e.g. iMac Low fc Klessenl (120041 ): lElmegreen &; Scalo 

(|2004j). Supernova are widely-believed to be the most importa nt source of turbulence for diffuse gas, 

and t his has been demonstrated by numerical simulations (e.g., Korpi et alll999 : de Avillez Breitschwerdt 
2005). In the outer galaxy where s tar formation rates are low, driving by the magnetorota tional 



instability may, however, dominate (jSellwood k, Balbus 



1999 



Piontek fc Ostrikeril2005l . 12007m . Spi- 
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ral shocks are al so effective in driving turbulence in the warm ISM, especially at high latitudes 
(|Kim et al.l 12008). The interaction between large-scale self-gravity, rotation , and shear can drive 
near-sonic turbulence at large scales (jWada et alj|2002l : iKim &: Ostrikerll2007l ). although the ampli- 
tude of this at scales less than the disk scale height may be modest. It is not known how effective 
these other mechanisms are for driving turbulence within GMCs, however, which are very dense 
and therefore present a small effective crossection to the diffuse ISM. 

3. Feedback and the Toomre Q parameter 

We have measured the Toomre parameter for each of our models, considering both the entire 
medium and just the dense gas, and comparing "thermal-only" with "turbulent +thermal" values. 
For all models, we find that the "turbulent +thermal" Q-value for dense gas is in the range 1 — 2, 
and is much greater than the thermal-only value. A further interesting point is that the turbulence 
level evidently adjusts with surface density in order to reach a marginally-stable state. In Series Q 
and R, which have f2/S = const., the velocity dispersions are relatively independent of S, yielding 
marginally-stable Q in the cold, dense gas. In Series K, which has constant £1 and therefore is highly 
unstable at large £ in the absence of turbulence, the velocity dispersions strongly increase with S, as 
a consequence of much higher levels of feedback activity. These higher velocity dispersions lift Q to 
near unity. Thus, our simulations give direct evidence of feedback lea ding to a self - regula ted quasi- 
steady state. We note that Q values vary for different components; IWada et al.l (J2002j) similarly 
found a large range of Q when measured in local patches within their disk simulations. 

Depending on what exactly is included in a model, the threshold fo r gravitational instabilit y 
in previous non-turbulent simulations is measured to be at Q ~ 1.5 (see lMcKee fc Ostrikerl 120071 ). 
which is similar to the values we find here when turbulence is included in Q. In the actively- 
star-forming regions o f galaxies, measured values of the Toomre parameter are not constant, but 
show a limited range (Marti n fc Kennicuttl 120011). Evidence for star formation thresho lds tied to 



Q are more mixed ( Martin Kennicuttl 



2001 



Boissier et al. 



2003; 



Boissier et al 



because star formation in outer disk s primarily takes place in spiral arms (jFerguson et al. 



2007) , possibly 



1998 



Thilker et al.l 120071 : iBush et al.ll2008l ) which strongly compress the gas above ambient densities. 



4. The virial ratio 

We measure the virial ratio 1Z (eq. 1551) for all our models, separating into different density 
regimes. We find that dense gas (n > 100 cm -3 ) generally has 1Z between 1 — 2, whereas lower- 
density gas has large values of 1Z. 1Z does not vary strongly with X in any of the series. In 
particular, we note that in spite of the large range of velocity dispersions in the dense gas in Series 
K at different S, 1Z varies only weakly with S. This indicates that feedback can effectively regulate 
the dynamics within massive, dense clouds, independent of the larger-s cale galactic environ ment. 
This is consistent with both older studies based on 12 CO obs ervations (jSolomon et al.lll987l ). and 
recent studies based on 13 CO observations ( Hever et al.l l2008). both of which find 1Z near 1 — 2 for 
Milky Way GMCs. Although masses based on CO are less cert ain in external gala xies, virial ratios 
also likely near unity for GMCs observed in the Local Group (IBolatto et al.l l2008) 
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5. Cloud surface densities 



We estimate the surface density of typical clouds by measuring the mass-weighted vertically- 
integrated column of gas; we define this as S c i 0U( j. For the values of the feedback threshold that we 
adopt, we find that £ c ioud is in the range 70 — 150 M pc -2 for most models. This is comparable 
to the ty pical GMC surface densities that are observed in the Mi l ky Way and in L ocal Group 
galaxies (jSolomon et al.lll987l : Isheth et al.1 12008 : iBolatto et al.1 12008 : iHever et aUbood ). Whether 
observational selection effects or physical processes impose a limited range of column densities 
for GMCs is an open question. For example, magnetic fields may impose a minimum surface 
density for formation o f gravitationally-bqund s tructures in the ISM, at a value E = B / (2ir^/G) = 
30 M pc 2 (S/10//G) (jMcKee fc Ostrikerl MXM ) ■ Here, we find that altering the volume heated 
in our HII region prescription does not appreciably change E c i ou d, but changing the gravitational 
potential threshold for star formation feedback does: when potential thresholds are low, S c i OU( j is 
also low. Comparison of cloud properties with observations may turn out to be a much more critical 
test of whether an ISM model is realistic than some other measures, such as the star formation 
rate. 

6. Dependence of EgpR on £ 

We obtain estimates of the dependence of surface star formation rate EsFR ° n g as surface 
density £ in our models by assuming that the timescale for star formation is proportional to the 
free-fall time in gas above some density threshold n t h = Pth/ A 4 - We compare results for two different 
threshold densities, nth = 100 cm -3 and re t h = 10 3 cm -3 . For Series Q and R (which most resemble 
real galaxies), we find (for either choice of threshold) relations that are well-described by power 
laws: Esfr oc £ 1+p with 1 +p = 1.2 — 1.4. In Series K, a power law is a less-good fit; the slope is 
also shallower (closer to unity). 

These resul ts are consistent with empirical Kennicutt-Schmidt laws ( Schmidt 19591. 19631; 



1989 



KennicuttJll998ah. which show similar values of 1+p when all the gas is includ ed in £ (e.g 



Kennicutt 



1998bl : IWong k Blitdl2002l : ISchuster et alJl2007l : lKennicutt et al.ll2007h . Re cent work has sug- 



gested that 1+ p is close to unity if just CO-emitting molecular gas is included (jBigiel et al. 2008 
in preparation! ) ; this implies all CO-emitting gas (most of which is at n = 10 2 — 10 3 cm -3 ) has 
the same star formation rate independent of galactic environment. Our prescription that the star 
formation rate per unit mass of dense gas is constant is equivalent to empirically finding 1 + p = 1 
if only molecular gas is included. 

It is encouraging that the results we find for Kennicutt-Schmidt relations in our disk models 
with feedback are compatible with observations. We also find, however, that star formation rates 
predicted from hydrostatic models are in fact similar to those predicted from the hydrodynamic 
models, with similar slopes at large E. This is true even though the hydrostatic models are not 
at all like real galaxies. Thus, one must be cautious in considering a numerical model successful 
if it yields reasonable star formation rates, since this can simply be a consequence of choosing 
reasonable initial conditions in a simulation. Indeed, a number of recent numerical studies have 
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found results similar to observec. 
they in cluded in the models (e.g 



20081 ). ISchave k Dalla Vecchia 



Kennicutt-Schmidt laws, regardles s of the detailed physics that 



Li et al.ll2006l : iTasker k Brvanll2006l . |200J; iRobertson k Kravtsov 



2009 ) have also recently emphasized that reproducing empirical 



star formation scaling laws is not by itself a critical test of an ISM model. 
7. Density- dependence of star formation efficiency 

The star formation efficiency per free-fall time can be defined locally as a function of threshold 
density by eg(pth) = tfr(/0th)SsFR./S(/) > p t h)', corresponding global measures can also be obtained. 
For a given true star formation rate, eg(p t h) oc tsf where the scaled star formation time tsf is 
shown for two different threshold densities in Fig. [TBI We find that tsf (or eg (/3th)) decreases 
with increasing p t h_- This is not a strong effect, however: it is less than a factor 2 for an order 
of magnitude difference in p t h_- In Series Q and R, the ratio of efficiencies at different threshold 
densities are also independent of £ (although this is not true for Series K). 



Krumholz k Tanl (|2007l ) recently compiled a range of observations of eg (which they refer to as 
SFRg). They point out that for threshold densities above ~ 100cm -3 , the value of eg does not vary 
strongly with density. They also find a smaller value for gas traced by HCN than for gas traced by 
CO, which since HCN has a higher critical density than CO is consistent with our finding that eg 
decreases with increasing p t h- 

Even though scaled star formation times at high density vary together independent of S, the 
free-fall time at the vertically-averaged density is not proportional to tsf- Instead, we find that 
TSF/te (psve) increases at increasing S. This implies that a prescription for star formation based 
on the mean density within one "average" scale height in a disk (or from a simulation that does 
not resolve high-density gas) would increasingly overestimate the star formation rate at high X. 
The same is true for the orbital time and the Jeans time: an assumption that the star formation 
rate varies oc £/i rb or oc S/tj would increasingly overestimate Ssfr at high £ compared to the 
value obtained from measuring the mass of gas in dense, gravitationally-bound regions. Thus, ISM 
models must resolve self-gravitating structures at scales less than the disk thickness in order to 
make accurate predictions of the star formation rate. 



We are grateful to the referee for a thorough reading and thoughtful set of comments that 
have helped us to clarify our presentation. This research was supported by grant NNG-05GG43G 
from NASA. Numerical computations were carried out on the OIT High Performance Computing 
Cluster (HPCC) and CTC cluster in the Department of Astronomy, at the University of Maryland. 
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Cartesian Disk Potential via Fourier Transforms 



In this section, we provide details for our method of solving Poisson's equation in Cartesian 
disk geometry. We generalize to the three-dimensional case. Thus, the solution is applicable to 
problems in which periodic boundary conditions in the horizontal (x-y) directions are assumed, 
and vacuum boundar y conditions a r e requ ired in the vertical (z) direction. Following Binney & 
Tremaine (1987) and iMiyama et al.1 (119871 ) . the gravitational potential at the location (x a ,Xb,x c ) 
on a regular Cartesian grid with dimensions (N x , N y , N z ) for integer indices (a, b, c) can be written 
as 



®(x a ,y b ,z c ) 



4ttG 

'N x Ny 



.v„_,A„- 1 .v : _, _ 2 ~i ( ^+ bn ^ 



E E E 

m=0 n=0 j=0 



N x Ny) ~ ( 7 .\n ( 7 7 .\ 



(Al) 



Here G m ^ n is the Green function of the Poisson equation for a horizontal sheet sinusoidal source 
charge and is written as 



G m ,n{z, Z ) — < 



1 



Oh 
1, 



exp( k m ^\z z |) (km,n 7^ 0), 



z — z 



{k"m,n — 0); 



(A2) 



where 



km,n 



2irn 



\TT) + V L 



1/2 



(A3) 



Note that up to a constant, the k m ^ n = solution is equal to the limit of the general solution for 
k m ,n 0. The coefficients p m>n (zj) are given by the discrete Fourier transforms of the density in 
the horizontal plane z = Zj, i.e., 



Pm,n( z j) = E E e 



v i v / am bn \ 

N x -iN y -i 2m 1 

— — N X Ny) 



vJ p(x a ,y b ,Zj). 



(A4) 



a=0 6=0 



Since z — z' in the Green function takes on values between — L z and L z , we can apply Fourier 
transforms on the domain (—L Z ,L Z ) to obtain 



1 

2N~ Z 



2N Z -1 

E ■ 



-2m 
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Substitution of equation (|A2p and evaluation of the sum shows that 



1 - exp (±i£tt - k m ^ n L z ) 

t"m,n,£= ; ^2 ■ v A 'J 



U2 -\- I ll 

If we now substitute into the original expression for <£, we obtain 



, ' am bn cl 
-2m — + — + 



47rG N x -lN y -l2N z -l g \N X N y 2N Z/ 

<t>(x a ,x b ,x c ) = - N Yl E E p L «:^ ( A8 ) 



where 



fcL^H^F) +(^) (A9) 
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and 



2N Z -1 27Ti jL 



/C^= E e 2 ^«fe)- (A10) 



Here, 



-corr = J Pm,nfe) for J = 0, N z - 1 

Pm ' n " \ -p m)n (^-ivje- fe -.^ for j = N z , 2N Z - 1 1 j 

That is, in addition to using the horizontal Fourier components of the original domain, we must 
define an image density in the augmented domain, based on the Pm,n values a distance L z away. 
After the corrected density is defined, a Fourier transform in z is taken in the usual way. 

The complete procedure for obtaining the gravitational potential is therefore as follows: First 
compute two-dimensional (x — y) Fourier components for each horizontal layer in z. Next define 
Pm V n( z j) anci take the Fourier transform in z to obtain p c ^ x n i- Finally, multiply by the Poisson 
kernel —AirG/k^ », and take the three-dimensional inverse transform. We note that even and odd 
terms in the sum on £ may also be combined, such that the forward and inverse transforms are 
both three dimensional; the Poisson kernel then is adapted to incorporate terms from the image 
charge. 

Since we use FFTs (via the Package FFTW) for all transforms the computational expense 
scales as 0(N log N), which is superior to the Green function method using direct summation. 
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Table 1: Dynamical Model Parameters 
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[a] Gas surface density, in Mq pc - 

[b] Orbital angular velocity, in km s _1 kpc -1 

[c] Orbital time, t or b = 2ir/Q, in units 10 8 yr 

[d] Duration of simulation, in units 10 s yr 

[e] Stellar density at midplane, in Mq pc~ 3 

[f] Vertical size of domain, in pc. Horizontal size is L x = 4L Z . 

[g] Average number density, n — H/(fj,L z ), in cm -3 

[h] Thermal conduction in units of Ko = 6.91 x 10 7 erg cm -1 s _1 K 

[i] Grid scale, in pc 

[j] Control parameter for photoheating region (see eq(34j 
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Table 2: Hydrostatic Model Parameters 



model 


S 


p* 


L z 


n 


K/Kq 


Az 


HSP8 • 


7.50 


0.035 


410 


0.525 


0.4 


0.40 


HSP11 • 


• 10.60 


0.07 


290 


1.05 


0.2 


0.28 


HSP15 • 


• 15.00 


0.14 


205 


2.1 


0.1 


0.20 


HSP21 • 


• 21.20 


0.28 


145 


4.2 


0.05 


0.14 


HSP42 • 


• 42.40 


1.12 


72.5 


16.8 


0.0125 


0.070 


HSC8 • 


7.50 


0.14 


353 


1.05 


0.1 


0.20 


HSC11 • 


• 10.60 






1.48 






HSC21 • 


• 21.20 






2.97 






HSC42 • 


• 42.40 






5.93 







Note. - 



All parameters and units are as in Table [T] 
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t = 329.5 Myr 
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X(pc) 
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Fig. 1. — A snapshot of Model Qll, with S = 10.6M Q pc" 2 : Upper panel shows temperature, and 
Lower panel shows density; both have logarithmic color scales as shown. Contours in the lower 
panel are loci of constant the relative potential, with solid and dashed curves for $W > and 
5>W < 0, respectively. The black asterisk at the potential minimum in the left part of the panel 
indicates the location of a "star zone" feedback center. The cold, dense gas is highly clumpy and 
filamentary, and is concentrated toward the midplane by self- and external gravity. 




Fig. 2. — A later-time (At= 38 Myr) snapshot of temperature (top) and density (bottom) from 
Model Qll, showing that individual structures have dynamically evolved. See Figure Q] legend for 
details. 
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t = 329.5 Myr t = 367.6 Myr 




log (n) log (n) 

Fig. 3. — Phase diagram for Model Qll snapshots shown in Figures [1] and [2 The gray scale shows 
the relative fraction of the total mass at each point in the phase plane. The solid line denotes the 
thermal equilibrium curve at the background heating rate (Gfuv = 1); the dashed curve shows 
thermal equilibrium in HII regions (Gruv = 10 3 )- The snapshot from Figure Q] is taken while an 
HII region is very young, and there is considerable gas on the hot, high-density branch; there is 
very little gas on the HII region branch from the second snapshot. Overall, gas is preferentially - 
but not exclusively - found near the thermal equilibrium curves. 
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Fig. 4. — Density distributions for Model Qll snapshot shown in Figure [TJ Solid and dashed lines 
denote mass- and volume- weighted density PDFs. The bimodal distribution of mass is character- 
istic of the modified two-phase warm/cold medium that develops for a bistable cooling function 
with a range of pressures (due to vertical stratification, self-gravity, and turbulent dynamics). 
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Fig. 5. — Time evolution of component energies in Model Qll: Solid, Long-dashed, and Short- 
dashed lines denote the thermal, kinetic, and potential energies, respectively, integrated over the 
domain. For the potential energy, we use the relative potential (see eqn. [23]) . 
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Fig. 6. — A snapshot of temperature (top) and density (bottom) in high surface density Model Q42. 
See Figure [H legend for details. 
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Fig. 7. — Mass fractions of gas components in all model series, as a function of surface density 
(Series Q, K, R) or angular velocity (Series S). Open symbols denote neutral gas at densities in the 
primarily-atomic regime; filled symbols denote neutral gas at densities in the primarily-molecular 
regime; asterisks denote photoheated gas corresponding to HII regions. 
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Fig. 8. — The mass-weighted surface density of gas (S c i ou d; see eq. H8|) as a function of area- 
weighted surface density E. Values of E c i ou d are similar to observed GMC surface densities, and 
relatively independent of environmental conditions. 
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Fig. 9. — Mean temperature in varying density ranges. Diffuse (WNM: open circles, CNM: open 
squares), dense (DM2: filled triangles, DM3: filled circles), and photoheated (HII: asterisks) com- 
ponents are defined as in Fig. [JJ 
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E[Solar Mass/pc 2 ] £2[km/s/kpc] 

Fig. 10. — Mean turbulent velocity dispersion in each gas component. Turbulence levels are rela- 
tively independent of £ in Series Q and R which have k and 17 oc S, but turbulence increases with 
£ in Series K, which has k and O constant for all models and is therefore highly unstable at large 
X. Turbulence levels decrease with increasing Q in series S, because angular momentum stabilizes 
against gravitational collapse and subsequent feedback. 
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Fig. 11. — Values of the Toomre Q parameter. We separately compute Q in four different ways: 
including all gas and either using the thermal (subscript "T") or thermal+turbulent velocity dis- 
persion to obtain Qt (total) and Q (total), respectively; and including just the (cold) high-density 
gas to obtain (thermal) Qt(^ > 100) and (thermal+turbulent) Q(n > 100). Feedback drives turbu- 
lence and regulates the thermal+turbulent Toomre parameter to reach levels near unity, especially 
for cold, dense gas. See text for details. 
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Fig. 12. — Virial ratio 1Z = 2E/\W\ (see text), measured separately for each gas component. High 
density gas (n > 100 cm" 3 ) is strongly gravitationally bound; low density gas is not. 
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Fig. 13. — Vertically-averaged density for all dynamic (open squares) and hydrostatic model series. 
Filled squares and circles denote hydrostatic series HSP and HSC, respectively. Turbulence reduces 
the mean density (averaged over the disk scale height) relative to the hydrostatic case. Mean density 
increases with S more rapidly in Series Q compared to Series R because the latter adopts the same 
stellar vertical gravity value for all models. 
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Fig. 14. — Ratio of mean free-fall to orbital time. Open squares denote hydro dynamic models, and 
filled squares denote the hydrostatic series (HSP and HSC). 
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Fig. 15. — Scaled star formation time, tsf = isF£ff(pth) = tff(pth)£/£(p > Pth), computed using 
two different density thresholds n t h = Pth/p- Results for dynamical models are shown with open 
triangles and circles; results for hydrostatic models are shown with filled triangles. For Series Q 
and R (which are most similar to real galaxies at varying radius), the two choices of threshold 
would yield consistent values of the star formation time tsp provided that the adopted efficiency 
decreases slightly with density, eff(n t h = 100) ~ (1.5 — 2)cg (n t h = 10 3 ). 
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Fig. 16. — Scaled star formation rates per unit area, EsFRAff(Pth) = S(p > Pth)Aff(pth)- Open 
triangles and circles are based on threshold density of nth = 10 2 cm" 3 and n t h = 10 3 cm" 3 , 
respectively. Filled triangles show the results for hydrostatic models. The wedge-shape marks 
indicate reference slopes of 1.0 and 1.5. Series Q and R show Kennicutt-Schmidt indices similar to 
observations. Efficiency parameters eg ~ 0.001 — 0.01 for n t h = 10 2 — 10 3 cm" 3 would be required 
to match the observed range of Esfr in this range of gas surface density E. 
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Fig. 17. — Ratio of scaled star formation time to orbital time, indicating that the ratio is not a 
constant independent of gas surface density S. 
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Fig. 18. — Ratio of scaled star formation time to Jeans time; the Jeans time is evidently not a good 
predictor of the star formation timescale. 
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Fig. 19. — Ratio of scaled star formation time to free fall time at large-scale mean density averaged 
over a disk scale height, ts(p avc ). An estimated star formation timescale proportional to the free-fall 
time at the vertically-averaged density (i.e. the "vertically-unresolved" limit) would increasingly 
underpredict the true star formation time at high S. 



